使用 R 計算邏輯回歸中的係數

在多元線性回歸中，可以使用以下公式找出係數。

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

例如：

> y <- c(9.3, 4.8, 8.9, 6.5, 4.2, 6.2, 7.4, 6, 7.6, 6.1)
> x0 <- c(1,1,1,1,1,1,1,1,1,1) 
> x1 <-  c(100,50,100,100,50,80,75,65,90,90)
> x2 <- c(4,3,4,2,2,2,3,4,3,2)
> Y <- as.matrix(y)
> X <- as.matrix(cbind(x0,x1,x2))

> beta = solve(t(X) %*% X) %*% (t(X) %*% Y);beta
        [,1]
x0 -0.8687015
x1  0.0611346
x2  0.9234254
> model <- lm(y~+x1+x2) ; model$coefficients
(Intercept)          x1          x2 
-0.8687015   0.0611346   0.9234254 

我想如何以相同的“手動”方式計算邏輯回歸的 beta。y 當然是 1 或 0。假設我使用的是帶有 logit 鏈接的二項式系列。

線性回歸模型中的 OLS 估計器很少有可以用封閉形式表示的性質，即不需要表示為函數的優化器。然而，它是一個函數的優化器——殘差平方和函數——並且可以這樣計算。

邏輯回歸模型中的 MLE 也是適當定義的對數似然函數的優化器，但由於它不能用於封閉形式的表達式，因此必須將其計算為優化器。

大多數統計估計器只能表示為適當構造的數據函數的優化器，稱為標準函數。這種優化器需要使用適當的數值優化算法。函數的優化器可以在 R 中使用optim()提供一些通用優化算法的函數或更專業的包之一（如optimx. 了解針對不同類型的模型和統計標準函數使用哪種優化算法是關鍵。

線性回歸殘差平方和

OLS 估計器被定義為著名的殘差平方和函數的優化器：

在像殘差平方和這樣的二次可微凸函數的情況下，大多數基於梯度的優化器都做得很好。在這種情況下，我將使用 BFGS 算法。

#================================================
# reading in the data & pre-processing
#================================================
urlSheatherData = "http://www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv"
dfSheather = as.data.frame(read.csv(urlSheatherData, header = TRUE))

# create the design matrices
vY = as.matrix(dfSheather['InMichelin'])
mX = as.matrix(dfSheather[c('Service','Decor', 'Food', 'Price')])

# add an intercept to the predictor variables
mX = cbind(1, mX)

# the number of variables and observations
iK = ncol(mX)
iN = nrow(mX)

#================================================
# compute the linear regression parameters as 
#   an optimal value
#================================================
# the residual sum of squares criterion function
fnRSS = function(vBeta, vY, mX) {
 return(sum((vY - mX %*% vBeta)^2))
}

# arbitrary starting values
vBeta0 = rep(0, ncol(mX))

# minimise the RSS function to get the parameter estimates
optimLinReg = optim(vBeta0, fnRSS,
                  mX = mX, vY = vY, method = 'BFGS', 
                  hessian=TRUE)

#================================================
# compare to the LM function
#================================================
linregSheather = lm(InMichelin ~ Service + Decor + Food + Price,
                   data = dfSheather)

這產生：

> print(cbind(coef(linregSheather), optimLinReg$par))
                   [,1]         [,2]
(Intercept) -1.492092490 -1.492093965
Service     -0.011176619 -0.011176583
Decor        0.044193000  0.044193023
Food         0.057733737  0.057733770
Price        0.001797941  0.001797934

邏輯回歸對數似然

Logistic回歸模型中MLE對應的準則函數是對數似然函數。

在哪裡是邏輯函數。參數估計是這個函數的優化器

我再次展示瞭如何optim()使用 BFGS 算法構建和優化標準函數。

#================================================
# compute the logistic regression parameters as 
#   an optimal value
#================================================
# define the logistic transformation
logit = function(mX, vBeta) {
 return(exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) )
}

# stable parametrisation of the log-likelihood function
# Note: The negative of the log-likelihood is being returned, since we will be
# /minimising/ the function.
logLikelihoodLogitStable = function(vBeta, mX, vY) {
 return(-sum(
   vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta)))
   + (1-vY)*(-log(1 + exp(mX %*% vBeta)))
   ) 
 ) 
}

# initial set of parameters
vBeta0 = c(10, -0.1, -0.3, 0.001, 0.01) # arbitrary starting parameters

# minimise the (negative) log-likelihood to get the logit fit
optimLogit = optim(vBeta0, logLikelihoodLogitStable,
                  mX = mX, vY = vY, method = 'BFGS', 
                  hessian=TRUE)

#================================================
# test against the implementation in R
# NOTE glm uses IRWLS: 
# http://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares
# rather than the BFGS algorithm that we have reported
#================================================
logitSheather = glm(InMichelin ~ Service + Decor + Food + Price,
                                 data = dfSheather, 
                        family = binomial, x = TRUE)

這產生

> print(cbind(coef(logitSheather), optimLogit$par))
                   [,1]         [,2]
(Intercept) -11.19745057 -11.19661798
Service      -0.19242411  -0.19249119
Decor         0.09997273   0.09992445
Food          0.40484706   0.40483753
Price         0.09171953   0.09175369

作為警告，請注意數值優化算法需要謹慎使用，否則您最終可能會遇到各種病態的解決方案。在您很好地理解它們之前，最好使用可用的打包選項，讓您專注於指定模型，而不是擔心如何數值計算估計值。

引用自：https://stats.stackexchange.com/questions/81000