[PYTHON] Implémenter un modèle de régression logistique en temps discret avec stan
TL;DR Si vous souhaitez effectuer une analyse du temps de survie où les covariables varient dans le temps, vous pouvez utiliser une technique appelée régression logistique en temps discret. Nous allons introduire la régression logistique en temps discret et l'implémenter avec stan.
Au début
Imaginez une situation où vous souhaitez estimer les facteurs d'obtention du diplôme sans redoubler un an après la transition des notes universitaires. Mettez les notes de l'élève $ i $ et la note $ t $ comme $ X_ {i, t} $, et la note t comme $ Y_ {i, t} $. Si vous pensez à un état dans lequel vous continuez à réussir une promotion comme une «survie», cela est saisi dans le cadre de l'analyse du temps de survie.
On s'attend à ce que l'ampleur des facteurs que $ X_ {i, t} $ aura sur la promotion variera considérablement en fonction de la note. Par exemple, il est difficile pour les élèves de première et de deuxième année de récupérer par la suite même s'ils n'obtiennent pas certains crédits, et les élèves de la troisième à la quatrième année doivent remplir les conditions d'obtention du diplôme, il est donc facile pour eux de redoubler l'année. Je vais.
La première chose qui vient à l'esprit lorsque l'on réfléchit au cadre permettant de savoir si l'on peut obtenir son diplôme sans redoubler d'un an est le modèle de risque proportionnel utilisé en pharmacie. Cependant, le modèle de hasard proportionnel suppose que les covariables ne changent pas avec le temps, donc si $ X_ {i, t} $ change avec le temps, vous devez l'examiner.
Par conséquent, nous utilisons un modèle de régression logistique en temps discret.
Modèle de régression logistique en temps discret
Dans l'analyse du temps de survie, la fonction de survie $ S (t) $ (probabilité de survivre pendant t période ou plus) et la fonction de risque $ h (t) $ (probabilité de survivre au point temporel t et de mourir au point temporel t) sont définies comme suit. Je vais.
Dans le modèle de régression logistique en temps discret, la fonction de risque est exprimée comme suit.
C'est la même chose qu'une régression logistique normale. $ \ Beta_t $ représente la différence d'efficacité des covariables en fonction du moment. Modifiez la forme de la fonction z (x, t) en fonction de votre situation.
Exemple de mise en œuvre de stan
Supposons que l'utilisateur continue de choisir l'un des choix multiples à tout moment. Estimez quelle option vous empêchera de partir le plus et comment l'effet de l'option changera à un moment donné.
modèle
stan_code = '''
//référence: 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; //Terme constant
}
model {
vector[ST] zeros = rep_vector(0.0, ST); //Vecteur à fixer à 0 lors de l'estimation du coefficient
vector[C] ones = rep_vector(1.0, C); //Matrice de sommation des colonnes
vector[ST] mu = alpha + (X_score .* append_col(zeros, beta[X_T, :])) * ones; //Vecteur de probabilité de retrait
for (st in 1:ST) {
target += bernoulli_logit_lpmf(Y[st] | mu[st]); //Modèle de classification binaire qui prédit le retrait
}
}
'''
Génération de données virtuelles
def logistic_func(array):
"""-∞~Réglez la valeur d'entrée de ∞ en fonction de la fonction logistique[0,1]Convertir et revenir"""
return 1/(1+np.exp(-array))
#Génération de données,Coupé à droite
S = 10000
C = 3
alpha = -1 #Taux de sortie par défaut avant multiplication par la fonction logistique, c'est environ 20%Sur
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 #Augmentez le nombre de colonnes pour correspondre à R,Ajuster l'efficacité du coefficient avec la dernière valeur
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() #Notez que cela stocke un tableau à l'intérieur
stan_dict["Y"] = list()
stan_dict["S"] = list() #Pour le débogage
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: #Jugement de survie
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: #Il sera interrompu lorsqu'il dépasse un certain niveau
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)
La vraie valeur de la version bêta a été définie comme ceci.
Le choix de l'option intermédiaire dans les premiers stades augmentera le taux de survie, mais dans les étapes finales, le choix de la bonne option augmentera le taux de survie par rapport à l'option intermédiaire.
La distribution du temps de survie ressemble à ceci.
Estimé
#Post-traitement des données transmises à stan
stan_dict["ST"] = np.sum(stan_dict["S"])
stan_dict["X_score"] = np.array(stan_dict["X_score"])
#Modéliser et déduire
sm = pystan.StanModel(model_code = stan_code)
fit = sm.sampling(stan_dict, iter=1000, chains=4)
print(fit)
Rhat ne dépasse pas 1,1, il semble donc converger.
Vérifiez le résultat
Vraie valeur de la bêta
estimation bêta
Les valeurs sont proches les unes des autres. Cependant, la précision de l'estimation est pire dans la seconde moitié de la période de survie.
référence
https://www.slideshare.net/akira_11/dt-logistic-regression
Recommended Posts