Jak symulować sztuczne dane dla regresji logistycznej?

45

Wiem, że brakuje mi czegoś w rozumieniu regresji logistycznej i naprawdę doceniłbym każdą pomoc.

O ile rozumiem, regresja logistyczna zakłada, że ​​prawdopodobieństwo wyniku „1” przy danych wejściowych jest liniową kombinacją danych wejściowych, przechodzącą przez funkcję odwrotnej logistyki. Jest to zilustrowane w następującym kodzie R:

#create data:
x1 = rnorm(1000)           # some continuous variables 
x2 = rnorm(1000)
z = 1 + 2*x1 + 3*x2        # linear combination with a bias
pr = 1/(1+exp(-z))         # pass through an inv-logit function
y = pr > 0.5               # take as '1' if probability > 0.5

#now feed it to glm:
df = data.frame(y=y,x1=x1,x2=x2)
glm =glm( y~x1+x2,data=df,family="binomial")

i pojawia się następujący komunikat o błędzie:

Komunikaty ostrzegawcze: 1: glm.fit: algorytm nie jest zbieżny 2: glm.fit: dopasowane prawdopodobieństwa wystąpiły 0 lub 1

Od jakiegoś czasu pracuję z R. wystarczy wiedzieć, że prawdopodobnie to ja jestem winien ... co się tutaj dzieje?

zorbar
źródło
2
Sposób, w jaki symulujesz swoje dane, wydaje mi się dziwny. Jeśli chcesz, dla alternatywnego bardziej standardowego sposobu, możesz zajrzeć tutaj: stats.stackexchange.com/questions/12857/…
ocram
@ocram: masz rację; to jest duplikat pytania!
user603,
2
Przeprowadziłem błędną symulację, jak wyjaśnił @ Stéphane Laurent. Problemem była jednak idealna separacja w regresji logistycznej , problem, którego nie znałem i raczej byłem zaskoczony, gdy się o nim dowiedziałem.
zorbar
@zorbar: to była moja odpowiedź na twoje pytanie (teraz usunięte).
user603
1
@ user603: Prawdopodobnie przegapiłem twoją odpowiedź; W każdym razie dzięki
zorbar,

Odpowiedzi:

63

yi1pr(i)

> set.seed(666)
> x1 = rnorm(1000)           # some continuous variables 
> x2 = rnorm(1000)
> z = 1 + 2*x1 + 3*x2        # linear combination with a bias
> pr = 1/(1+exp(-z))         # pass through an inv-logit function
> y = rbinom(1000,1,pr)      # bernoulli response variable
> 
> #now feed it to glm:
> df = data.frame(y=y,x1=x1,x2=x2)
> glm( y~x1+x2,data=df,family="binomial")

Call:  glm(formula = y ~ x1 + x2, family = "binomial", data = df)

Coefficients:
(Intercept)           x1           x2  
     0.9915       2.2731       3.1853  

Degrees of Freedom: 999 Total (i.e. Null);  997 Residual
Null Deviance:      1355 
Residual Deviance: 582.9        AIC: 588.9 
Stéphane Laurent
źródło
Masz rację - przegapiłem ten krok. bardzo dziękuję za Twoją pomoc!
zorbar,
1
Miałem pytanie dotyczące sposobu symulacji danych. Symulując dane dla regresji liniowej, symulujemy również szum (\ epsilon). Rozumiem, że funkcja logistyczna jest funkcją oczekiwania, która sama eliminuje hałas. Czy to jest powód, dla którego nie masz hałasu w swoim Z?
Sam
5
mean response+noise
@ Dokładnie StéphaneLaurent. Całkowicie rozumiem. Dziękuję bardzo za odpowiedź.
Sam
2

Regresja logistyczna nadaje się do dopasowania, jeśli jako wartości docelowe podano prawdopodobieństwa lub proporcje, a nie tylko wyniki 0/1.

import numpy as np
import pandas as pd
def logistic(x, b, noise=None):
    L = x.T.dot(b)
    if noise is not None:
        L = L+noise
    return 1/(1+np.exp(-L))

x = np.arange(-10., 10, 0.05)
bias = np.ones(len(x))
X = np.vstack([x,bias]) # Add intercept
B =  [1., 1.] # Sigmoid params for X

# True mean
p = logistic(X, B)
# Noisy mean
pnoisy = logistic(X, B, noise=np.random.normal(loc=0., scale=1., size=len(x)))
# dichotomize pnoisy -- sample 0/1 with probability pnoisy
dichot = np.random.binomial(1., pnoisy)

pd.Series(p, index=x).plot(style='-')
pd.Series(pnoisy, index=x).plot(style='.')
pd.Series(dichot, index=x).plot(style='.')

Tutaj mamy trzy potencjalne cele regresji logistycznej. pktóra jest wartością rzeczywistą / docelową / prawdopodobieństwem, pnoisyktóra jest równa p przy normalnym szumie dodanym do logarytmicznej skali szans, i dichotktóra jest nieprzyzwoita traktowana jako parametr dwumianowego pliku PDF i na tej podstawie pobierana jest próbka. Powinieneś przetestować wszystkie 3 - zauważyłem, że niektóre implementacje LR typu open source nie pasują p.

W zależności od aplikacji możesz preferować cierpliwość.

W praktyce powinieneś również rozważyć, w jaki sposób hałas może się kształtować w docelowej aplikacji i spróbować to naśladować.

użytkownik48956
źródło