Hitung koefisien dalam regresi logistik dengan R

18

Dalam regresi linier berganda dimungkinkan untuk mengetahui koefisien dengan rumus berikut.

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

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

Contohnya:

> 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

Saya ingin cara menghitung dalam cara "manual" yang sama beta untuk regresi logistik. Di mana tentu saja y akan menjadi 1 atau 0. Dengan asumsi saya menggunakan keluarga binomial dengan tautan logit.

logistic regression-coefficients S12000
sumber

1

Pertanyaan yang sebenarnya Anda tanyakan telah diajukan di stats.stackexchange.com/questions/949/… . Pertanyaan yang sepertinya ingin Anda tanyakan adalah jawabannya di sini.

Whuber

26

Estimator OLS dalam model regresi linier sangat jarang memiliki properti yang dapat direpresentasikan dalam bentuk tertutup, yang tanpa perlu dinyatakan sebagai pengoptimal fungsi. Namun demikian, ini adalah pengoptimal dari suatu fungsi - jumlah sisa dari fungsi kuadrat - dan dapat dihitung dengan demikian.

MLE dalam model regresi logistik juga merupakan pengoptimal dari fungsi log-likelihood yang sesuai, tetapi karena itu tidak tersedia dalam ekspresi bentuk tertutup, itu harus dihitung sebagai pengoptimal.

Sebagian besar penaksir statistik hanya dapat diekspresikan sebagai pengoptimal fungsi yang dibangun dengan tepat dari data yang disebut fungsi kriteria. Pengoptimal semacam itu membutuhkan penggunaan algoritma pengoptimalan numerik yang tepat. Pengoptimal fungsi dapat dihitung dalam R menggunakan optim()fungsi yang menyediakan beberapa algoritma optimasi tujuan umum, atau salah satu paket yang lebih khusus seperti optimx. Mengetahui algoritma pengoptimalan mana yang akan digunakan untuk berbagai jenis model dan fungsi kriteria statistik adalah kuncinya.

Regresi linear jumlah residual kuadrat

Estimator OLS didefinisikan sebagai pengoptimal dari fungsi residu jumlah kuadrat yang terkenal:

\begin{aligned} \hat{β} & = \arg min_{β} {(Y - X β)}^{'} (Y - X β) \\ = (X^{'} X)^{- 1} X^{'} Y \end{aligned}

$\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}$

Dalam kasus fungsi cembung yang dapat dibedakan dua kali seperti jumlah kuadrat residu, sebagian besar pengoptimal berbasis gradien melakukan pekerjaan dengan baik. Dalam hal ini, saya akan menggunakan algoritma 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)

Ini menghasilkan:

> 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

Logistik regresi log-kemungkinan

Fungsi kriteria yang sesuai dengan MLE dalam model regresi logistik adalah fungsi log-likelihood.

\begin{aligned} \log L_{n} (β) & = \sum_{i = 1}^{n} (Y_{i} \log Λ (X_{i}^{'} β) + (1 - Y_{i}) \log (1 - Λ (X_{i}^{'} β))) \end{aligned}

$\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}$ mana adalah fungsi logistik. Estimasi parameter adalah pengoptimal dari fungsi ini

Λ (k) = 1 / (1 + \exp (- k))

$\Lambda(k) = 1/(1+ \exp(-k))$

\hat{β} = \arg max_{β} \log L_{n} (β)

$\hat{\boldsymbol{\beta}} = \arg\max_{\boldsymbol{\beta}}\log L_n(\boldsymbol{\beta})$

Saya menunjukkan bagaimana membangun dan mengoptimalkan fungsi kriteria menggunakan optim()fungsi sekali lagi menggunakan algoritma 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)

Ini menghasilkan

> 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

Sebagai peringatan, perhatikan bahwa algoritma pengoptimalan numerik memerlukan penggunaan yang cermat atau Anda dapat berakhir dengan segala macam solusi patologis. Sampai Anda memahami mereka dengan baik, yang terbaik adalah menggunakan opsi paket yang tersedia yang memungkinkan Anda untuk berkonsentrasi menentukan model daripada khawatir tentang bagaimana menghitung estimasi secara numerik.

tchakravarty
sumber

1

Kerja bagus @tchakravarty, fungsi log-likelihood dapat disederhanakan dengan menggunakan-sum(vY%*%(mX%*%vBeta)-log(1+exp(mX%*%vBeta)))

Mamoun Benghezal

11

Anda tidak bisa sampai di sana dari sini. Solusi untuk model linear umum dan model logistik muncul dari penyelesaian persamaan kemungkinan maksimum masing-masing, tetapi hanya model linier yang memiliki solusi bentuk tertutup.

Jika Anda membaca buku McCullagh dan Nelder, Anda dapat mempelajari bagaimana solusinya diperoleh dalam kasus logistik (atau model umum lainnya). Akibatnya, solusi dihasilkan secara iteratif, di mana setiap iterasi melibatkan penyelesaian regresi linier tertimbang. Bobot sebagian tergantung pada fungsi tautan.

Placidia
sumber

2

atau cari di web untuk "kuadrat terkecil berulang berulang" ...

Ben Bolker

Atau langsung di sini: stats.stackexchange.com/questions/236676/…

kjetil b halvorsen

Hitung koefisien dalam regresi logistik dengan R

Jawaban:

Regresi linear jumlah residual kuadrat

Logistik regresi log-kemungkinan