[PYTHON] Implementieren Sie mit stan ein zeitdiskretes logistisches Regressionsmodell
TL;DR Wenn Sie eine Überlebenszeitanalyse mit Kovariaten durchführen möchten, die sich im Laufe der Zeit ändern, können Sie eine Technik verwenden, die als zeitdiskrete logistische Regression bezeichnet wird. Wir werden eine zeitdiskrete logistische Regression einführen und mit stan implementieren.
Am Anfang
Stellen Sie sich eine Situation vor, in der Sie die Faktoren für den Abschluss abschätzen möchten, ohne ein Jahr nach dem Übergang der Hochschulnoten zu wiederholen. Setzen Sie die Noten des Schülers $ i $ und die Note $ t $ als $ X_ {i, t} $ und die Note t als $ Y_ {i, t} $. Wenn Sie sich einen Zustand, in dem Sie weiterhin erfolgreich befördert werden, als "Überleben" vorstellen, wird dies im Rahmen der Überlebenszeitanalyse erfasst.
Es wird erwartet, dass die Größe der Faktoren, die $ X_ {i, t} $ bei der Beförderung haben wird, je nach Besoldungsgruppe stark variieren wird. Zum Beispiel ist es für Erst- und Zweitklässler schwierig, sich danach zu erholen, selbst wenn sie keine Credits erhalten, und 3. bis 4. Klasse müssen die Abschlussanforderungen erfüllen, so dass es für sie einfach ist, das Jahr zu wiederholen. Ich werde.
Das erste, was Ihnen einfällt, wenn Sie über den Rahmen nachdenken, ob Sie Ihren Abschluss machen können, ohne ein Jahr zu wiederholen, ist das in der Pharmazie verwendete Proportional-Hazard-Modell. Das Proportional-Hazard-Modell geht jedoch davon aus, dass sich die Kovariaten im Laufe der Zeit nicht ändern. Wenn sich also $ X_ {i, t} $ im Laufe der Zeit ändert, müssen Sie dies überprüfen.
Daher verwenden wir ein zeitdiskretes logistisches Regressionsmodell.
Zeitdiskretes logistisches Regressionsmodell
In der Überlebenszeitanalyse werden die Überlebensfunktion $ S (t) $ (Überlebenswahrscheinlichkeit für t Zeitraum oder mehr) und die Gefahrenfunktion $ h (t) $ (Überlebenswahrscheinlichkeit bis zum Zeitpunkt t und Sterben zum Zeitpunkt t) wie folgt definiert. Ich werde.
Im zeitdiskreten logistischen Regressionsmodell wird die Gefährdungsfunktion wie folgt ausgedrückt.
Dies entspricht einer normalen logistischen Regression. $ \ Beta_t $ repräsentiert den Unterschied in der Wirksamkeit von Kovariaten in Abhängigkeit vom Zeitpunkt. Ändern Sie die Form der Funktion z (x, t) entsprechend Ihrer Situation.
Implementierungsbeispiel von stan
Angenommen, der Benutzer wählt zu einem bestimmten Zeitpunkt immer wieder eine von mehreren Möglichkeiten. Schätzen Sie, welche Option Sie daran hindert, die Option am meisten zu verlassen, und wie sich die Wirkung der Option zu einem bestimmten Zeitpunkt ändert.
Modell-
stan_code = '''
//Referenz: https://www.slideshare.net/akira_11/dt-logistic-regression
data {
int T_max;
int ST;
int C;
int X_T[ST];
matrix[ST, C] X_score;
int Y[ST];
}
parameters {
matrix[T_max, C-1] beta;
real alpha; //Konstante Laufzeit
}
model {
vector[ST] zeros = rep_vector(0.0, ST); //Vektor, der bei der Schätzung des Koeffizienten auf 0 festgelegt werden soll
vector[C] ones = rep_vector(1.0, C); //Matrix zum Summieren von Spalten
vector[ST] mu = alpha + (X_score .* append_col(zeros, beta[X_T, :])) * ones; //Entnahmewahrscheinlichkeitsvektor
for (st in 1:ST) {
target += bernoulli_logit_lpmf(Y[st] | mu[st]); //Binäres Klassifizierungsmodell, das den Entzug vorhersagt
}
}
'''
Virtuelle Datengenerierung
def logistic_func(array):
"""-∞~Stellen Sie den Eingabewert von ∞ entsprechend der Logistikfunktion ein[0,1]Konvertieren in und zurück"""
return 1/(1+np.exp(-array))
#Datengenerierung,Rechts abschneiden
S = 10000
C = 3
alpha = -1 #Die Standard-Exit-Rate vor dem Multiplizieren mit der Logistikfunktion beträgt ca. 20%Über
T_max = 6
beta = np.array([[0,i,j] for i, j in zip(list(range(-3, 3)), list(range(-3, 3))[::-1])]) * 0.2 #Erhöhen Sie die Anzahl der Spalten entsprechend R.,Passen Sie die Effektivität des Koeffizienten an den letzten Wert an
stan_dict = dict()
stan_dict["T_max"] = T_max
stan_dict["C"] = C
stan_dict["class"] = list()
stan_dict["rate"] = list()
stan_dict["X_T"] = list()
stan_dict["X_score"] = list() #Beachten Sie, dass hierdurch ein Array gespeichert wird
stan_dict["Y"] = list()
stan_dict["S"] = list() #Zum Debuggen
for s in range(S):
idx = 0
class_ = np.random.choice(list(range(C)), size=1)[0]
x_score = np.zeros((C))
x_score[score] = 1.0
rate = logistic_func(alpha+beta[idx,score])
stan_dict["class"].append(score)
stan_dict["rate"].append(rate)
stan_dict["X_T"].append(idx+1)
stan_dict["X_score"].append(x_score)
while True: #Überlebensurteil
if int(np.random.binomial(n=1, p=rate, size=1)):
y = 1
stan_dict["Y"].append(y)
stan_dict["S"].append(idx+1)
break
elif idx >= T_max-1: #Es wird eingestellt, wenn es ein bestimmtes Niveau überschreitet
y = 0
stan_dict["Y"].append(y)
stan_dict["S"].append(idx+1)
break
y = 0
stan_dict["Y"].append(y)
idx += 1
score = np.random.choice(list(range(R)), size=1)[0]
x_score = np.zeros((R))
x_score[score] = 1.0
rate = logistic_func(alpha+beta[idx,score])
stan_dict["class"].append(score)
stan_dict["rate"].append(rate)
stan_dict["X_T"].append(idx+1)
stan_dict["X_score"].append(x_score)
Der wahre Wert von Beta wurde so eingestellt.
Die Auswahl der mittleren Option in den frühen Stadien erhöht die Überlebensrate, aber in der letzten Phase erhöht die Auswahl der richtigen Option die Überlebensrate als die mittlere Option.
Die Verteilung der Überlebenszeit sieht so aus.
Schätzen
#Nachbearbeitung von Daten an Stan
stan_dict["ST"] = np.sum(stan_dict["S"])
stan_dict["X_score"] = np.array(stan_dict["X_score"])
#Modellieren und schließen
sm = pystan.StanModel(model_code = stan_code)
fit = sm.sampling(stan_dict, iter=1000, chains=4)
print(fit)
Rhat überschreitet nicht 1,1, also scheint es zu konvergieren.
Überprüfen Sie das Ergebnis
Wahrer Wert von Beta
Beta-Schätzung
Die Werte liegen nahe beieinander. Die Schätzgenauigkeit ist jedoch in der zweiten Hälfte der Überlebenszeit schlechter.
Referenz
https://www.slideshare.net/akira_11/dt-logistic-regression
Recommended Posts