Расчет коэффициентов регрессии logisti c вручную - PullRequest
Новая
       100

Расчет коэффициентов регрессии logisti c вручную

2 голосов
/ 17 июня 2020

Я пытаюсь оценить логистическую c регрессию в R, рассчитывая все вручную. Я могу создать функцию lo git и loglikelihood, однако я не могу решить ее с помощью некоторого нелинейного решателя

Я хотел бы спросить совета 

df <- read_csv("http://courses.atlas.illinois.edu/spring2016/STAT/STAT200/RProgramming/data/Default.csv")
df


df$default = ifelse(df$default == "Yes", 1, 0)

logit <- function(x, b0, b1) {
  1/(1 + exp(-b0 - b1*x))
}



Loglikel <- function(y, x, b0, b1) {
  b0 = rep(b0, length(y))
  b1 = rep(b1, length(y))
  p <- logit(x, b0, b1)
  sum(y*log(p)  + (1 - y)*log(1-  p))
}


Loglikel(df$default, df$balance, -10, 0.005)


library(stats4)

mle(Loglikel, 
    start = list(b0 = 0, b1 = 0), 
    fixed = list(y = df$default, x = df$balance))

1 Ответ

1 голос
/ 17 июня 2020

Я взял ваш код и немного изменил его, чтобы передавать параметры в виде вектора: 

df <- read_csv("http://courses.atlas.illinois.edu/spring2016/STAT/STAT200/RProgramming/data/Default.csv")
df$default <- ifelse(df$student == "Yes", 1, 0)

logit <- function(x, b0, b1) {
  1/(1 + exp(-b0-b1*x))
}
Loglikel <- function(par, y, x){
  p <- logit(x, par[1], par[2])
  sum(y*log(p)  + (1-y)*log(1-p))
}

Теперь мы готовы использовать нелинейный решатель (например, nlm) для получения оценок для параметров: 

nlm_fit <- nlm(Loglikel, p = c(-2,0.001), x=df$balance, y=df$default)

, что дает 

> nlm_fit
...
$estimate
[1] -2.0002960 -0.2666521
...

nlm, использует решатель типа Ньютона-Рафсона для минимизации MLE. В то же время glm использует алгоритм наименьших квадратов с итеративным взвешиванием, что означает, что выходные данные glm и nlm не должны согласовываться: 

glm_fit <- glm(default ~ balance, family = binomial(link="logit"), data = df)

> glm_fit

Call:  glm(formula = default ~ balance, family = binomial(link = "logit"), 
    data = df)

Coefficients:
(Intercept)      balance  
 -1.7004224    0.0009409

Проверьте эту ссылку , это дает хорошее представление о том, что происходит под капотом glm.

...