TL;DR

공변수가 시변의 생존 시간 분석을 하려면 이산 시간 논리 회귀 방법을 사용할 수 있다.
이산 시간 물류 회귀를 소개하고stan을 통해 실현한다.

개시하다

대학 성적의 추이를 통해 유급되지 않고 졸업하는 주요 원인을 추측하다.
학생$i$, 학년$t달러의 성적$X{i,t}$,학년 t가 $Y로 순조롭게 승급할 수 있을까{i,t} 저기에 놓으세요.
진급이 지속적으로 성공한 상태를'생존'으로 이해하면 생존 시간 분석의 틀을 통해 포착할 수 있다.
$X_{i,t}달러가 진급에 주는 주요 원인의 크기는 학년에 따라 큰 차이가 있을 수 있습니다.예컨대 1·2학년은 학점을 다소 못 받더라도 나중에 만회할 수 있어 유급이 어렵고, 3·4학년은 졸업 요건을 충족해야 해 유급이 용이하다.
유급 없이 졸업할 수 있는 틀을 고려할 때 가장 먼저 떠오르는 것은 약학 등이 사용하는 비례위험모델이다.
그러나 비례의 위험 모델은 공통 변수가 시간에 따라 변하지 않는다고 가정하기 때문에 $X{i,t}달러는 시간의 변화에 따라 꼼꼼하게 고려해야 한다.
따라서 이산 시간 논리 회귀 모델을 사용한다.

이산 시간 논리 회귀 모델

생존 시간 분석에서 생존 함수 $S(t)$(t 기간 이상 생존 확률)와 위험 함수 $h(t)$(t)$(t)$(시간 t까지 생존할 때 사망할 확률)는 다음과 같이 정의된다.
$$
S(t) = P(T>t) =\int_{x}^{\infty} f(t)dt\\
h(t) =\lim_{\Delta\rightarrow\infty}\frac{P(t\leq T < t+\Delta t | T\geq t)}{\Delta t}\\
$$
이산 시간 논리 회귀 모델에서 위험 함수는 다음과 같다.
$$
h(x,t) =\frac{1}{1 +\exp(-z(x,t))}\\
z(x,t) =\alpha +\beta_t * X_{t}\\
$$
이것은 일반적인 물류 회귀와 같다.
$\beta_시점에 따른 공통 변수 효과의 차이를 나타냅니다.
상황에 따라 z(x, t) 함수의 모양을 바꾸십시오.

stan의 실현 예시

사용자가 항상 여러 옵션 중 하나를 선택한다고 가정합니다.
어떤 선택항이 가장 이탈을 방지할 수 있는지, 선택항의 효과에 따라 시간대에 어떻게 변화하는지 추측한다.

모델

stan_code = '''
// 参考: 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; // 定数項
}

model {
  vector[ST] zeros = rep_vector(0.0, ST); // 係数推定時に、0に固定するためのベクトル
  vector[C] ones = rep_vector(1.0, C); // 列に対して和をとるための行列
  vector[ST] mu = alpha + (X_score .* append_col(zeros, beta[X_T, :])) * ones; // 離脱確率ベクトル

  for (st in 1:ST) {
    target += bernoulli_logit_lpmf(Y[st] | mu[st]); // 離脱を予測する2値分類モデル
  }
}
'''

가상 데이터 생성

def logistic_func(array):
    """-∞~∞の入力値をlogistic functionに則って[0,1]に変換して返す"""
    return 1/(1+np.exp(-array))

# データの生成, 右打ち切りにする
S = 10000
C = 3
alpha = -1 # logistic関数掛ける前のデフォルトの離脱率、これで約20%くらい
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 # Rに合わせて列を増やすこと, 係数の効き具合は最後の値で調整
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() # これは中にarrayを格納するので注意
stan_dict["Y"] = list()
stan_dict["S"] = list() # デバッグ用

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: # 生存判定
        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: # 一定以上になったら打ち切り
            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)

베타의 실제 값은 이렇게 설정되어 있다.

중간의 옵션을 초기에 선택하면 생존율이 높아지지만 마지막에 오른쪽의 옵션을 선택하면 중간의 옵션보다 생존율이 높아진다.
생존 기간의 분포는 이렇다.

어림잡다

# stanに与えるデータの後処理
stan_dict["ST"] = np.sum(stan_dict["S"])
stan_dict["X_score"] = np.array(stan_dict["X_score"])

# モデル立てて推論
sm = pystan.StanModel(model_code = stan_code)
fit = sm.sampling(stan_dict, iter=1000, chains=4)
print(fit)

Rhat이 1.1을 넘지 않았기 때문에 수렴될 수 있다.

결과를 확인하다

beta의 실제 값

beta의 추정치

가격 차이가 많지 않다.그러나 생존 기간 후반부의 추정 정밀도는 떨어진다.

참고 자료