полиномиальная логит-регрессия вручную в R

Я пытаюсь реализовать полиномиальную регрессию (пакет mlogit или multinom) в R с кодами и оптимизацией (без использования пакетов).

rm(list= ls())
data = read.table("~/Desktop/R Code/textfiles/keane.csv", sep = ",",header = T)

data1 = data[,c("educ","exper", "expersq", "black", "status")]
data1 = na.omit(data1)

data2 = as.matrix(data1)


y_1 = rep(0, nrow(data2))
y_2 = rep(0, nrow(data2))
y_3 = rep(0, nrow(data2))


data2 = cbind(data2[,1:5], y_1, y_2, y_3)

data2[,6] = ifelse(data2[,5] == 1, 1, 0)
data2[,7] = ifelse(data2[,5] == 2, 1, 0)
data2[,8] = ifelse(data2[,5] == 3, 1, 0)


int = rep(1, nrow(data2)) #intercept

data2 = cbind(int, data2[,c(1:4,6,7,8)]) 


X = as.matrix(data2[, c(1:5)])
y_1 = as.matrix(data2[, 6]) #replace y values(status = 1)
y_2 = as.matrix(data2[, 7]) #replace y values(status = 2)
y_3 = as.matrix(data2[, 8]) #replace y values(status = 3)


Y = cbind(y_1, y_2, y_3) 


##beta


beta = solve(t(X) %*% X) %*% t(X) %*% Y #LPM coefficient 


logit.nll = function (beta, X, Y) {
    
    
    
    
    P = as.matrix(rowSums(exp(X %*% beta))); #Sum_(h=1)^3 exp(X * Beta_(h))
    
    
    
    Pr_1 = exp(X %*% beta[,2])/(1 + P); #P(y = 2 | X)
    Pr_2 = exp(X %*% beta[,3])/(1 + P); #P(y = 3 | X)
    
    Pr_0 = 1/(1+P);#P(y = 1 | X)
    
 
 
 (colSums(Y[,1] * log(Pr_0)) + colSums(Y[,2] * log(Pr_1)) + colSums(Y[,3] * log(Pr_2))) #log-likelihood
    
    
    
}

optim(beta, logit.nll, X = X, Y = Y, method = "BFGS")

когда я делаю этот код, он дает мне сообщение об ошибке в X %*% beta: несоответствующие аргументы. Мой подход может быть в корне неправильным, или реализация функции логарифмического правдоподобия неверна. Могу ли я получить помощь, чтобы исправить этот код?


r logistic-regression multinomial
person John    schedule 01.09.2020    source источник
comment
Без доступа к keane.csv мы не можем воспроизвести проблему, и это затрудняет диагностику. Можете ли вы разместить эти данные где-нибудь? (Обратите внимание, что вы можете использовать dput() в R, а затем скопировать/вставить сюда вывод dput() как один из способов предоставления данных).   -  person DanY    schedule 01.09.2020

Ответы (1)


arrow_upward
1
arrow_downward

Не очень хорошо знаком с вашей оптимизацией svm или с тем, что вы пытаетесь сделать, ваша ошибка связана с optim работой с вектором. Вам нужно преобразовать его в матрицу внутри функции, допустим, ваши данные такие:

set.seed(111)
data = iris

X = model.matrix(~.,data=data[,1:4])
Y = model.matrix(~0+Species,data=data)

beta = solve(t(X) %*% X) %*% t(X) %*% Y

Теперь мы добавляем матричную часть, также обратите внимание, что по умолчанию оптим выполняет минимизацию (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html), поэтому вам нужно вернуть отрицательное значение логарифмической вероятности:

logit.nll = function (beta, X, Y) {
    
    beta = matrix(beta,ncol=3)    
    P = as.matrix(rowSums(exp(X %*% beta))); #Sum_(h=1)^3 exp(X * Beta_(h))    
    Pr_1 = exp(X %*% beta[,2])/(1 + P); #P(y = 2 | X)
    Pr_2 = exp(X %*% beta[,3])/(1 + P); #P(y = 3 | X)    
    Pr_0 = 1/(1+P);#P(y = 1 | X)
     
LL = (colSums(Y[,1] * log(Pr_0)) + colSums(Y[,2] * log(Pr_1)) + colSums(Y[,3] * log(Pr_2))) #log-likelihood
print(LL)
return(-LL)

}

res = optim(beta, logit.nll, X = X, Y = Y, method = "BFGS")

res
$par
             Speciessetosa Speciesversicolor Speciesvirginica
(Intercept)      -2.085162         15.040679        -27.60634
Sepal.Length     -4.649971         -8.971237        -11.43702
Sepal.Width      -9.286757         -5.016616        -11.69764
Petal.Length     12.803070         17.125483         26.55641
Petal.Width       6.025760          3.342659         21.63200
person StupidWolf    schedule 01.09.2020
comment
Большое спасибо, просто добавление бета-термина в функцию работает отлично! - person John; 01.09.2020
comment
Почему это дает такие разные результаты, чем nnet::multinom(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data = iris)? - person invictus; 15.01.2021