Я пытаюсь воспроизвести результаты функции multinom () с функцией optim () в R, но это не дает тех же результатов. Что случилось?

Question

Я хочу воспроизвести результаты функции multinom () с функцией optim () в R, но она не дает таких же результатов. Что было не так?

Сначала я импортировал данные publi c как «ml».

require(foreign)
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")

Коды для получения сводной статистики данных «ml» и результатов ниже:

with(ml, table(ses,prog))
with(ml, do.call(rbind,tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))

        prog
ses      general academic vocation
  low         16       19       12
  middle      20       44       31
  high         9       42        7

                M       SD
general  51.33333 9.397775
academic 56.25714 7.943343
vocation 46.76000 9.318754

Коды для получения результатов из функции multinom (), которая выполняет полиномиальную логистическую c регрессию, и следующие результаты ниже:

library(nnet)
ml$prog2 <- relevel(ml$prog, ref = "academic")
ml_pckg <- multinom(prog2 ~ write + ses, data = ml)
summary(ml_pckg)

Call:
multinom(formula = prog2 ~ write + ses, data = ml)

Coefficients:
         (Intercept)      write  sesmiddle    seshigh
general     2.852198 -0.0579287 -0.5332810 -1.1628226
vocation    5.218260 -0.1136037  0.2913859 -0.9826649

Std. Errors:
         (Intercept)      write sesmiddle   seshigh
general     1.166441 0.02141097 0.4437323 0.5142196
vocation    1.163552 0.02221996 0.4763739 0.5955665

Residual Deviance: 359.9635 
AIC: 375.9635 

Код для получения статистики z и результатов ниже:

z <- summary(ml_pckg)$coefficients/summary(ml_pckg)$standard.errors
z

         (Intercept)     write  sesmiddle   seshigh
general     2.445214 -2.705562 -1.2018081 -2.261334
vocation    4.484769 -5.112689  0.6116747 -1.649967

Затем я написал код для воспроизведения результатов, приведенных выше. Я сгенерировал фиктивные переменные для категориальных зависимых / независимых переменных, как показано ниже:

ml$prog_academic <- ifelse(ml$prog == "academic", 1, 0)
ml$prog_general <- ifelse(ml$prog == "general", 1, 0)
ml$prog_vocational <- ifelse(ml$prog == "vocational", 1, 0)

ml$ses_low <- ifelse(ml$ses == "low", 1, 0)
ml$ses_middle <- ifelse(ml$ses == "middle", 1, 0)
ml$ses_high <- ifelse(ml$ses == "high", 1, 0)

Я сгенерировал один вектор для умножения с перехватом и подмножеством записи, ses_middle и ses_high для объясняющей переменной. ses_low здесь является базовым. Я назначил эти ковариаты новому фрейму данных с именем «X».

one <-as.data.frame(rep(1,200)) 
covar <- ml[,c(7,19,20)]
X <- data.frame(one,covar) #200*4

Затем я создал другой фрейм данных для зависимых переменных с именем «Y», который состоит из prog_general и prog_vocational. Здесь prog_academi c - это базовая линия.

Y <- ml[,16:17] #200*2

Я установил начальное значение параметров, аналогичное результатам функции mlo git (), чтобы функция оптимизации сходилась.

B_0 <- c(3, -0.1, -0.5, -1, 5, -0.1, 0.2, -1) #8*1 #initial value as vector

Здесь я обращаюсь к документу , чтобы найти вероятность полиномиальной логистической c регрессии. Вероятность находится в уравнении 31 на странице 12. Я обнаружил, что вторая часть уравнения также должна быть суммирована по i.

Я создал пустую матрицу «xb», чтобы включить в него часть

xb <- matrix(0, nrow=200, ncol=2) #200*2

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

mlogit <- function(B){
  B <- matrix(B, nrow=2, ncol=4, byrow=T) 
  for (i in 1:nrow(xb)){  #i is the dimension of individual: 200
    for (j in 1:ncol(xb)){  #j is the dimension of dependant variables -1 (categorical): 2
      xb[i,j] <- sum(X[i,]*B[j,]) #200*2
    }
  }
  
  exp <- exp(xb) #200*2
  sumexp <- rowSums(exp) #200*1
  sumexp <- as.numeric(sumexp)
  
  yxb <- Y*xb #200*2
  sumyxb <- sum(yxb)
  
  ll <-  sumyxb-sum(log(1+sumexp))
  -ll
}

mlogit_result <- optim(par = B_0, fn = mlogit)
mlogit_result

Результаты ниже:

$par
[1]  0.05325004 -0.01417267 -0.64375499 -0.96137147  6.33471560 -0.86154161  0.92387035 -0.65728823

$value
[1] 103.7692

$counts
function gradient 
     353       NA 

$convergence
[1] 0

$message
NULL

Если результаты соответствуют результатам функции mlo git (), $ par должен быть таким, как показано ниже:

2.852198 -0.0579287 -0.5332810 -1.1628226  5.218260 -0.1136037  0.2913859 -0.9826649

Я снова и снова просматривал свой код и функцию правдоподобия, но не нашел здесь ничего неправильного. . Я думаю, что, возможно, начальный параметр установлен неправильно или у созданной мной функции возникла проблема.

Может ли кто-нибудь дать мне какие-нибудь предложения по решению этой проблемы?

Спасибо!