Oblicz współczynniki w regresji logistycznej z R

W wielokrotnej regresji liniowej można znaleźć współczynnik za pomocą następującego wzoru.

$b = (X'X)^{-1}(X')Y$

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

Na przykład:

> 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 

Chciałbym obliczyć w ten sam „ręczny” sposób beta dla regresji logistycznej. Gdzie oczywiście y byłoby 1 lub 0. Zakładając, że używam rodziny dwumianowej z łączem logit.

Estymator OLS w modelu regresji liniowej jest dość rzadki, gdy ma właściwość, że może być reprezentowany w postaci zamkniętej, to znaczy bez potrzeby wyrażania go jako optymalizator funkcji. Jest to jednak optymalizator funkcji - funkcja resztkowej sumy kwadratów - i może być obliczony jako taki.

MLE w modelu regresji logistycznej jest również optymalizatorem odpowiednio zdefiniowanej funkcji prawdopodobieństwa logarytmicznego, ale ponieważ nie jest dostępny w postaci wyrażenia zamkniętego, należy go obliczyć jako optymalizator.

Większość estymatorów statystycznych można wyrazić jedynie jako optymalizatory odpowiednio skonstruowanych funkcji danych zwanych funkcjami kryterialnymi. Takie optymalizatory wymagają zastosowania odpowiednich algorytmów optymalizacji numerycznej. Optymalizatory funkcji można obliczyć w R za pomocą optim()funkcji, która udostępnia niektóre algorytmy optymalizacji ogólnego przeznaczenia lub jeden z bardziej specjalistycznych pakietów, takich jak optimx. Kluczowa jest znajomość algorytmu optymalizacji, który należy zastosować dla różnych typów modeli i funkcji kryteriów statystycznych.

Regresja liniowa resztkowa suma kwadratów

Estymator OLS jest definiowany jako optymalizator dobrze znanej resztkowej funkcji kwadratów:

$$\begin{align} \hat{\boldsymbol{\beta}} &= \arg\min_{\boldsymbol{\beta}}\left(\boldsymbol{Y} - \mathbf{X}\boldsymbol{\beta}\right)'\left(\boldsymbol{Y} - \mathbf{X}\boldsymbol{\beta}\right) \\ &= (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\boldsymbol{Y} \end{align}$$

W przypadku podwójnie różniczkowej funkcji wypukłej, takiej jak resztkowa suma kwadratów, większość optymalizatorów opartych na gradiencie dobrze sobie radzi. W takim przypadku będę używać algorytmu 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)

Daje to:

> 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

Logarytmiczna regresja prawdopodobieństwa

Funkcja kryterium odpowiadająca MLE w modelu regresji logistycznej jest funkcją logarytmu prawdopodobieństwa.

$$\begin{align} \log L_n(\boldsymbol{\beta}) &= \sum_{i=1}^n \left(Y_i \log \Lambda(\boldsymbol{X}_i'\boldsymbol{\beta}) + (1-Y_i)\log(1 - \Lambda(\boldsymbol{X}_i'\boldsymbol{\beta}))\right) \end{align}$$ gdzie to funkcja logistyczna. Szacunki parametrów są optymalizatorami tej funkcji $\Lambda(k) = 1/(1+ \exp(-k))$ $$\hat{\boldsymbol{\beta}} = \arg\max_{\boldsymbol{\beta}}\log L_n(\boldsymbol{\beta})$$

optim()Pokazuję, jak skonstruować i zoptymalizować funkcję kryterialną, używając funkcji jeszcze raz wykorzystującej algorytm 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)

To daje

> 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

Ostrzegamy, że algorytmy optymalizacji numerycznej wymagają starannego użycia lub możesz skończyć z różnego rodzaju rozwiązaniami patologicznymi. Dopóki nie zrozumiesz ich dobrze, najlepiej skorzystać z dostępnych opcji pakietowych, które pozwalają skoncentrować się na określeniu modelu, zamiast martwić się o to, jak obliczyć liczbowo szacunki.

Nie możesz się stąd dostać. Rozwiązania zarówno ogólnego modelu liniowego, jak i modelu logistycznego wynikają z rozwiązania odpowiednich równań maksymalnego prawdopodobieństwa, ale tylko model liniowy ma rozwiązanie w postaci zamkniętej.

Jeśli zajrzysz do książki McCullagha i Neldera, możesz dowiedzieć się, jak uzyskać rozwiązania w przypadku logistycznym (lub innym uogólnionym modelu). W efekcie rozwiązania są tworzone iteracyjnie, przy czym każda iteracja wymaga rozwiązania ważonej regresji liniowej. Wagi zależą częściowo od funkcji łącza.