Model logistyczny
MŁ
Przykładowe użycie podstawowych funkcji i objaśnienie, co jest czym
Najpierw prześledzimy działanie funkcji glm z podstawowego pakietu stats na prostym przykładzie, w którym jest tylko jedna zmienna objaśniająca, potem na trochę trudniejszym (więcej iksów).
Przykład prosty
Losuję przykładowe dane (nie jest bardzo istotne jak, ale jakby ktoś był ciekawy, to losuję z modelu logistycznego, żeby potem wszystko ładnie dało się na tych danych pokazać):
set.seed(123)
X <- rnorm(100)
Y <- as.factor(as.numeric(exp(3 + 5*X)/(1 + exp(3 + 5*X)) < runif(100, 0, 1)))
df_prosty <- data.frame(y=Y, x=X)Ważne jest to, że \(y\) jest jako factor. To może się zmieniać w zależności od pakietu, z którego korzystamy.
head(df_prosty)## y x
## 1 0 -0.56047565
## 2 1 -0.23017749
## 3 0 1.55870831
## 4 0 0.07050839
## 5 0 0.12928774
## 6 0 1.71506499str(df_prosty)## 'data.frame': 100 obs. of 2 variables:
## $ y: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 2 1 1 ...
## $ x: num -0.5605 -0.2302 1.5587 0.0705 0.1293 ...Dopasowuję model:
mod <- glm(y~x, data=df_prosty, family=binomial(link="logit"))Model logistyczny jest uogólnionym modelem liniowym (generalized linear model) i ta część binomial(link = "logit") do tego się odnosi. Parametr family: binomial mówi nam o rozkładzie zmiennej objaśnianej - u nas \(y\) ma rozkład \(0-1\), zaś link: logit mówi nam o funkcji łączącej, tzn. po jakim przekształceniu zakładamy liniową zależność.
Po lewej stronie logitowa funkcja łącząca, po prawej kombinacja liniowa: \[\log \left( \frac{\hat{p}(1|x)}{1 - \hat{p}(1|x)} \right) = \beta_0 + \beta' x.\]
\(\hat{p}(1|x)\) oznacza estymowane prawdopodobieństwo tego, że obserwacja należy do klasy \(1\) pod warunkiem, że wartość zmiennej objaśniającej dla tej obserwacji to \(x\) (bardziej formalnie można to zapisać tak: \(\hat{P}(Y_i = 1|X_i = x)\), gdzie \(i\) to indeks obserwacji).
Tutaj definiuję funkcję logit i odwrotną do niej:
logit <- function(x){
log(x / (1 - x))
}
# funkcja odwrotna do logit:
invlogit <- function(x){
exp(x) / (1 + exp(x))
}
# sprawdzenie:
logit(invlogit(0.4))## [1] 0.4A to jest dopasowany model. Poniżej wytłumaczenie, co jest czym :) Poniżej wypisuję wszystko na różne sposoby, żebyście widzieli, co się z czym łączy.
summary(mod)##
## Call:
## glm(formula = y ~ x, family = binomial(link = "logit"), data = df_prosty)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.73823 -0.22082 -0.02910 0.02325 1.92270
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9808 0.7708 -3.867 0.00011 ***
## x -5.9816 1.5280 -3.915 9.05e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.652 on 99 degrees of freedom
## Residual deviance: 41.742 on 98 degrees of freedom
## AIC: 45.742
##
## Number of Fisher Scoring iterations: 8Wyestymowany wektor \((\beta_0, \beta)\)
mod$coefficients ## (Intercept) x
## -2.980846 -5.981568Tak otrzymujemy estymowane prawdopodobieństwa \(\hat{p}(1|x)\) (wypisuję dla pięciu pierwszych obserwacji):
mod$fitted.values[1:5]## 1 2 3 4 5
## 5.918642e-01 1.674200e-01 4.531896e-06 3.221441e-02 2.288351e-02invlogit(mod$coefficients[1] + mod$coefficients[2]*df_prosty$x[1:5])## [1] 5.918642e-01 1.674200e-01 4.531896e-06 3.221441e-02 2.288351e-02# jesli nie podajemy argumentu newdata,
# to predict jest wykonywany dla zbioru,
# ktorego uzywalismy przy dopasowywaniu modelu
predict(mod, type="response")[1:5] ## 1 2 3 4 5
## 5.918642e-01 1.674200e-01 4.531896e-06 3.221441e-02 2.288351e-02prawdopodobienstwo1 <- predict(mod, type="response")
prawdopodobienstwo0 <- 1 - prawdopodobienstwo1W funkcji predict, jeśli ustawimytype="link", to zwracane będą wartości kombinacji liniowej (\(\beta_0 + \beta' x\))…
mod$coefficients[1] + mod$coefficients[2]*df_prosty$x[1:5]## [1] 0.3716772 -1.6040236 -12.3043656 -3.4025967 -3.7541893predict(mod, type="link")[1:5]## 1 2 3 4 5
## 0.3716772 -1.6040236 -12.3043656 -3.4025967 -3.7541893…czyli inaczej logarytm ilorazu szans tzn. \(\log\frac{\hat{p}(1|x)}{1 - \hat{p}(1|x)}\)
log(prawdopodobienstwo1/prawdopodobienstwo0)[1:5]## 1 2 3 4 5
## 0.3716772 -1.6040236 -12.3043656 -3.4025967 -3.7541893Ostatni możliwy type w predict.glm robi to:
predict(mod, type="terms")[1:5]## [1] 3.8932922 1.9175914 -8.7827505 0.1190184 -0.2325743mod$coefficients[2]*(df_prosty$x[1:5] - mean(df_prosty$x))## [1] 3.8932922 1.9175914 -8.7827505 0.1190184 -0.2325743A tak dopasowanie modelu logistycznego wygląda na wykresie (punkty z danych są na czarno, wyestymowane pstwa na różowo, żeby widać było lepiej wyróżniłam jeden punkt - czerwono-niebieski):
y <- as.numeric(as.character(df_prosty$y)) # pozbywam sie factora
x <- df_prosty$x
plot(y ~ x, pch=16)
points(x, prawdopodobienstwo1, pch=16, col="magenta", cex=0.5)
points(x[1], y[1], pch=16, col="red", cex=1.5)
points(x[1], prawdopodobienstwo1[1], pch=16, col="blue", cex=1.5)
abline(v=x[1], lty=2)
A tutaj na różowo i na niebiesko zaznaczona jest funkcja \(f(x) = \hat{p}(1|x) = \frac{\exp(\beta_0 + \beta x)}{1+\exp(\beta_0 + \beta x)}\):
y <- as.numeric(as.character(df_prosty$y)) # pozbywam sie factora
x <- df_prosty$x
plot(y ~ x, pch=16)
x.grid <- seq(min(x), max(x), length.out = 100)
prawdopodobienstwa.grid <- predict(mod, newdata=data.frame(x=x.grid), type="response")
lines(x.grid, prawdopodobienstwa.grid, col="magenta", lwd=2)
f <- function(x){
exp(mod$coefficients[1] + mod$coefficients[2]*x)/(1+exp(mod$coefficients[1] + mod$coefficients[2]*x))
}
lines(x.grid, f(x.grid), col="blue", lty=2)
Żeby dostać klasyfikator, musimy wybrać punkt odcięcia, tzn. wartość prawdopodobieństwa, od której będziemy klasyfikować obserwacje do klasy 1 (a pozostałe do 0). Może to być np. \(1/2\):
klasa <- ifelse(prawdopodobienstwo1 > 0.5, 1, 0)
tabela_reklasyfikacji <- table(prawdziwaKlasa=y, predykcja=klasa)
tabela_reklasyfikacji## predykcja
## prawdziwaKlasa 0 1
## 0 68 5
## 1 7 20accuracy <- sum(diag(tabela_reklasyfikacji)) / sum(tabela_reklasyfikacji) # to wszystko na zbiorze treningowym - przydalby się testowy
accuracy## [1] 0.88Do porównywania modeli zagnieżdżonych używa się statystyki testowej \[\Lambda = -2\log \frac{L(\omega)}{L(\Omega)} = 2 (l(\omega) - l(\Omega)),\] gdzie \(L(\omega)\) jest wiarogodnością (czyli wartość funkcji wiarogodności policzonej w punkcie odpowiadającym estymatorom największej wiarogodności w modelu \(\omega\)) modelu mniejszego, \(L(\Omega)\) wiarogodnością modelu większego, zaś \(l\) to log-wiarogodność. Ważna w modelu jest dewiancja, którą definiuje się w następujący sposób: \[dev_{\textrm{model}} = 2(l(\textrm{model wysycony}) - l(\textrm{model})).\] Intuicyjnie im większa dewiancja dla modelu, tym model gorzej dopasowany, bo bardziej różni się od modelu wysyconego. U nas \(l(\textrm{model wysycony}) = 0\), więc:
mod$deviance## [1] 41.74197-2*logLik(mod)## 'log Lik.' 41.74197 (df=2)Przyjrzyjmy się dokładniej funkcji log-wiarogodności w modelu logistycznym. Wzór jest taki:
\[L(\beta|y,x) = \prod_{i=1}^{n} \hat{p}(1|x_i)^{y_i}(1-\hat{p}(1|x_i))^{1-y_i},\] gdzie \(\hat{p}(1|x_i) = \exp(\beta_0 + \beta'x)/(1+\exp(\beta_0 + \beta'x))\), zatem \[l(\beta|y,x) = \sum_{i=1}^{n} \left(y_i\log(\hat{p}(1|x_i)) + (1-y_i)\log(1 - \hat{p}(1|x_i))\right).\]
Po przełożeniu na kod, wygląda to tak:
loglikelihood <- function(y, x, beta){
y <- as.numeric(as.character(y))
beta_0 <- beta[1]
beta <- beta[-1]
prawdopodobienstwo1 <- invlogit(beta_0 + t(beta)%*%x)
sum(y*log(prawdopodobienstwo1) + (1 - y)*log(1 - prawdopodobienstwo1))
}
loglikelihood(y=df_prosty$y, x=df_prosty$x, beta=mod$coefficients)## [1] -20.87098logLik(mod)## 'log Lik.' -20.87098 (df=2)Podobną do \(R^2\) miarą dopasowania modelu jest: \[1 - \frac{dev_{model}}{dev_{model0}},\] gdzie \(dev_{model0}\) jest dewiancją modelu pustego.
mod_pusty <- glm(y ~ 1, data=df_prosty, family=binomial(link="logit"))
1 - mod$deviance/mod_pusty$deviance## [1] 0.6421661 - mod$deviance/mod$null.deviance## [1] 0.642166Do porównania modeli zagnieżdżonych potrzebujemy trudniejszego przykładu. Na prostszym przykładzie omówimy jeszcze test Walda. z value to statystyka testowa w tym teście.
summary(mod)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.980846 0.7707673 -3.867374 1.100134e-04
## x -5.981568 1.5279940 -3.914654 9.053387e-05Powstaje w ten sposób:
summary(mod)$coefficients[2, 1] / summary(mod)$coefficients[2, 2]## [1] -3.914654summary(mod)$coefficients[2, 3]## [1] -3.914654zaś p-wartość liczona jest tak:
2*pnorm(summary(mod)$coefficients[2, 3])## [1] 9.053387e-05summary(mod)$coefficients[2, 4]## [1] 9.053387e-05Trudniejszy przykład (modele zagnieżdżone, wybór zmiennych)
Zmienna objaśniana zależy tylko od zmiennych \(x1-x5\), zmienne \(x6-x10\) są w modelu zbędne.
set.seed(123)
X <- matrix(rnorm(10000), 1000, 10)
colnames(X) <- paste0("x", 1:10)
term <- exp(3 + c(2, 1, -1, 0.5, -0.5)%*%t(X[,1:5]))
Y <- as.factor(as.numeric(term/(1 + term) < runif(100, 0, 1)))
df_trudniejszy <- data.frame(y=Y, X)head(df_trudniejszy)## y x1 x2 x3 x4 x5 x6
## 1 1 -0.56047565 -0.99579872 -0.5116037 -0.15030748 0.1965498 -0.4941739
## 2 0 -0.23017749 -1.03995504 0.2369379 -0.32775713 0.6501132 1.1275935
## 3 0 1.55870831 -0.01798024 -0.5415892 -1.44816529 0.6710042 -1.1469495
## 4 0 0.07050839 -0.13217513 1.2192276 -0.69728458 -1.2841578 1.4810186
## 5 0 0.12928774 -2.54934277 0.1741359 2.59849023 -2.0261096 0.9161912
## 6 0 1.71506499 1.04057346 -0.6152683 -0.03741501 2.2053261 0.3351310
## x7 x8 x9 x10
## 1 -0.6992281 -1.6180367 0.5110004 1.9315759
## 2 0.9964515 0.3791812 1.8079928 -0.6164747
## 3 -0.6927454 1.9022505 -1.7026150 -0.5625675
## 4 -0.1034830 0.6018743 0.2874488 -0.9899631
## 5 0.6038661 1.7323497 -0.2691142 2.7312276
## 6 -0.6080450 -0.1468696 -0.3795247 -0.7216662str(df_trudniejszy)## 'data.frame': 1000 obs. of 11 variables:
## $ y : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 1 ...
## $ x1 : num -0.5605 -0.2302 1.5587 0.0705 0.1293 ...
## $ x2 : num -0.996 -1.04 -0.018 -0.132 -2.549 ...
## $ x3 : num -0.512 0.237 -0.542 1.219 0.174 ...
## $ x4 : num -0.15 -0.328 -1.448 -0.697 2.598 ...
## $ x5 : num 0.197 0.65 0.671 -1.284 -2.026 ...
## $ x6 : num -0.494 1.128 -1.147 1.481 0.916 ...
## $ x7 : num -0.699 0.996 -0.693 -0.103 0.604 ...
## $ x8 : num -1.618 0.379 1.902 0.602 1.732 ...
## $ x9 : num 0.511 1.808 -1.703 0.287 -0.269 ...
## $ x10: num 1.932 -0.616 -0.563 -0.99 2.731 ...mod_duzy <- glm(y ~ ., data=df_trudniejszy, family=binomial(logit))
mod_maly <- glm(y ~ x1 + x2 + x3 + x4 + x5, data=df_trudniejszy, family=binomial(logit))
mod_bardzo_maly <- glm(y ~ x1 + x2, data=df_trudniejszy, family=binomial(logit))Test: \[ H_0: \textrm{model } \omega \textrm{ jest adekwatny} \\ H_1: \textrm{model } \Omega \textrm{ jest adekwatny} \] Model \(\omega\) ma \(q\) parametrów, model \(\Omega\) ma \(p > q\) parametrów. Statystyka testowa: \[D = dev_{\omega} - dev_{\Omega} \sim \chi_{p-q} \textrm{ przy } n \to \infty.\]
# test 1
stat1 <- mod_maly$deviance - mod_duzy$deviance
p <- 10
q <- 5
pchisq(stat1, p - q, lower.tail = FALSE)## [1] 0.3554668# test 2
q_bardzo_maly <- 2
stat2 <- mod_bardzo_maly$deviance - mod_maly$deviance
pchisq(stat2, q - q_bardzo_maly, lower.tail = FALSE)## [1] 3.808829e-14W teście 1 porównujemy model pełny (zmienne \(x1 - x10\)) z modelem ze zmiennymi \(x1 - x5\). P-wartość wychodzi duża - nie odrzucamy hipotezy zerowej o adekwatności modelu mniejszego. W teście 2 porównujemy model ze zmiennymi \(x1 - x5\) z modelem ze zmienymi \(x1 - x2\). Odrzucamy hipotezę zerową, tzn. co najmniej jedna ze zmiennych \(x3 - x5\) jest istotna.
To samo inaczej:
anova(mod_maly, mod_duzy, test="Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ x1 + x2 + x3 + x4 + x5
## Model 2: y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 994 607.53
## 2 989 602.00 5 5.5227 0.3555anova(mod_bardzo_maly, mod_maly, test="Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ x1 + x2
## Model 2: y ~ x1 + x2 + x3 + x4 + x5
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 997 673.09
## 2 994 607.53 3 65.559 3.809e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1To teraz selekcja zmiennych. Wszystko dzieje się bardzo podobnie, jak w przypadku funkcji lm. Pokażę na przykładzie AIC. \[AIC = dev_{\omega} + 2(p+1)\]
AIC(mod_duzy)## [1] 624.0035mod_duzy$deviance + 2*(p + 1)## [1] 624.0035# Nie chce wypisywac wszystkiego, dlatego trace=0
step(mod_duzy, trace=0)##
## Call: glm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x10, family = binomial(logit),
## data = df_trudniejszy)
##
## Coefficients:
## (Intercept) x1 x2 x3 x4 x5
## -2.4709 -1.6701 -0.7374 0.5586 -0.4778 0.5169
## x10
## 0.1829
##
## Degrees of Freedom: 999 Total (i.e. Null); 993 Residual
## Null Deviance: 911.8
## Residual Deviance: 604.5 AIC: 618.5Inne przydatne pakiety
Pakiet biglm
Do dużych danych (raczej tak dużych danych, żeby to się przydało, na zajęciach nie będzie, ale może warto wiedzieć, że jest)
library(biglm)
df_prosty$y <- as.numeric(as.character(df_prosty$y))
mod.biglm <- biglm::bigglm(y~x, data=df_prosty, family=binomial(logit), maxit=50)
summary(mod.biglm)## Large data regression model: bigglm(y ~ x, data = df_prosty, family = binomial(logit), maxit = 50)
## Sample size = 100
## Coef (95% CI) SE p
## (Intercept) -2.9808 -4.5224 -1.4393 0.7708 1e-04
## x -5.9816 -9.0376 -2.9256 1.5280 1e-04Pakiet glmnet
Możemy chcieć dopasować współczynniki modelu logistycznego korzystając z penalizowanej funkcji wiarogodności (czyli z karą).
library(glmnet)
x <- as.matrix(df_trudniejszy[,2:11])
y <- df_trudniejszy$y
# bez kary
mod.glmnet <- glmnet(x, y, family="binomial", lambda=0)
# z karą dla siatki lambd
mod.glmnet.lambda <- glmnet(x, y, family="binomial")Pakiet mlr
Pakiet, który ujednolica wywoływanie m.in. różnych klasyfikatorów i ma składnię podobną do pythonowego scikit-learna.
df_prosty$y <- as.factor(df_prosty$y)
library(mlr)
dane_treningowe <- makeClassifTask(data = df_prosty, target = "y")
jaki_model <- makeLearner("classif.logreg", predict.type = "response")
mod.mlr <- train(jaki_model, dane_treningowe)
predykcja <- predict(mod.mlr, dane_treningowe)
tabela_reklasyfikacji.mlr <- table(as.data.frame(predykcja)[,2:3])
tabela_reklasyfikacji.mlr## response
## truth 0 1
## 0 68 5
## 1 7 20