4.10 Else

4.10.1 Hw4. Rasch Model

$$ \[\begin{align} L(\theta, \beta) &= \prod_{k=1}^n \prod_{i=1}^p \left\{ \dfrac{\exp(\theta_k + \beta_i)}{1 + \exp(\theta_k + \beta_i)}\right\}^{y_{ki}} \left\{ \dfrac{1}{1 + \exp(\theta_k + \beta_i)}\right\}^{1-y_{ki}} \\ \pi (\theta, \beta \vert y) &= \pi(\theta) \pi(\beta) \ast \prod_{k=1}^n \prod_{i=1}^p \left\{ \dfrac{\exp(\theta_k + \beta_i)}{1 + \exp(\theta_k + \beta_i)}\right\}^{y_{ki}} \left\{ \dfrac{1}{1 + \exp(\theta_k + \beta_i)}\right\}^{1-y_{ki}} \end{align}\] $$

\[ 0<\left\{ \dfrac{\exp(\theta_k + \beta_i)}{1 + \exp(\theta_k + \beta_i)}\right\}^{y_{ki}} <1, \; \; \; \; \; 0<\left\{ \dfrac{1}{1 + \exp(\theta_k + \beta_i)}\right\}^{1-y_{ki}}<1 \]

underflow problem. log 취하면 해결.

\[ \pi(\theta_k) \sim N(0, \sigma^2), \; \; \; \; \; \pi(\sigma^2 ) \sim IG(0.001, 0.001) \]

update $\theta_k , k=1, \cdots, n$”

\[ \log (r) = \log \pi (\theta_k ' \vert y, \beta^{(t)}. \theta_{-k}^{(t)} - \log \pi (\theta^{(t)} \vert y, \beta^{(t)}. \theta_{-k}^{(t)} \]

if $\log U < min(\log(r), 0))$, accept. else, reject.

4.10.2 DA) Example: MVN

for DA 알고리즘, I-step과 P-step이 존재.

4.10.2.1 1. I-Step

$$ \[\begin{alignat}{4} Y_2 \vert Y_1 &\sim N && \Big( \mu_2 + \Sigma_{21} \Sigma_{11}^{-1} (Y_1 - \mu_1) , && \Sigma_{22} - \Sigma_{21}\Sigma_{11}^{-1} \Sigma_{12} \Big) \\ Y_{i, mis} \vert Y_{i, obs}, \mu, \Sigma &\sim N_{dim(Y_{mis}^{(i)})} && \Big( \mu_{mis}^{(i)} + \Sigma_{mis, obs}^{(i)} \Sigma_{obs, obs}^{-1} (Y_{i, obs} - \mu_{i, mis}^{(i)}) , && \Sigma_{mis,mis}^{(i)} - \Sigma_{mis,obs}^{(i)}[\Sigma_{obs,obs}^{(i)}]^{-1} \Sigma_{obs,mis}^{(i)} \Big) \end{alignat}\] $$

상기의 conditional pdf로 우리는 $Y_{i, mis}$를 impute 가능.

for $i=1, \cdots, n$, $Y_{i, mis} \vert Y_{i, obs}$ 에서 $Y_{i, mis}$를 draw.

4.10.2.2 2. P-Step

베이지안 분석을 위해선 prior가 필요. 여기서 prior는 이하로 설정하자. $q$는 known integer이며, $q=p$인 상황에 이는 $\Sigma$에 대한 Jefferey’s prior.

\[ \pi (\mu, \Sigma) \propto \vert \Sigma \vert^{-\tfrac{q+1}{2}} \]

위와 같이 식들을 구성하였을 때, com 데이터에 대한 posterior distribution $\pi(\mu, \Sigma \vert Y_1 , \cdots, Y_n)$는 이하와 같이 characterized 가능. 이는 inverse-Wishart 분포.

\[ \Sigma \vert Y_1 , \cdots, Y_n \sim \dfrac{1}{\vert \Sigma \vert^{-\tfrac{q+n}{2}}} \exp \left\{ -\dfrac{1}{2} tr\left( \Sigma^{-1} S \right) \right\} \]

이렇게 획득해온 패러미터들을 사용해 $\mu$의 post를 구하면 이하와 같다.

\[ \mu \vert \Sigma, Y_1 , \cdots, Y_n \sim N_p \left( \bar Y \dfrac {1}{n} \Sigma \right) \]

$\Sigma \vert Y_1 , \cdots, Y_n$에서 $\Sigma$ 를 생산

이후 $\mu \vert \Sigma, Y_1 , \cdots, Y_n$에서 $\mu$ 를 생산

4.10.3 Bayesian adaptive clinical trial with delayed outcomes

4.10.3.1 Continual Reassessment Method

Clinical Trial: Toxicity -> Efficacy -> Confirmation

희귀병 케이스에서는 도즈 레벨을 1~n까지 정해둔 후, 샘플을 slice 하여 1번 subsample에 도즈1 투입. 유효하면 (1차시에 3명 투입했다고 하고 그중에 1명이 독성 나왔으면 독성 확률은 1/3. 해당 여부로 도즈2로 넘어갈 것인지를 판단) 2투입. 2에서 문제 생기면 1로 복귀하고 2번 subsample에 도즈1 투입해봐서 유효한지 검증. 이렇게 모든 서브샘플에 도즈레벨 오가면서 투입해보고 최적 도즈레벨 결정.

이때 CRM을 시작하기 전 대략적으로 이정도의 도즈레벨이 최적 도즈레벨일 것이라는 예측 (Skeleton)을 정하고 CRM을 시작함.

delay outcome

이전 환자들의 evaluation이 끝나기 전에 (evaluation period가 경과하기 전이나, 결과가 나오기 전에) 환자 풀이 증가하는 상황

이 상황에서는 관측이 더 된 환자보다 덜 된 환자에서 outcome이 발생할 확률이 높음. 9개월 누워있던 놈보다 2개월 누워있던 놈이 12개월 경과 전에 뭔가 변화를 보이기 쉽다는 소리.

이때 아직 결과를 관찰하지 못한 환자들을 mis로 지정. 이 상황은 누워있던 기간이 결과 발생 여부라는 variant와 직결되어 있으므로 NMAR. 그러니까, 여기서 결과값은 outcome이 발생했는지 안했는지, 그리고 variable은 환자나 누워있던 기간.

위에서 언급했듯 누워있던 기간이 짧으면 변이확률 높음. 따라서 각각에 대해 다른 survival function을 적용하여 각각의 다른 확률 뽑아낸 후 이거 기반으로 DA 진행하면 해결.

4.10.4 NMAR의 종류

$m_i$는 missing indicator. $Y_i$가 mis면 1.

\[ f(M, Y \vert X, \theta, \psi) = \prod_{i=1}^n f(m_i , y_i \vert x_i , \theta, psi) \]

interested in direct relationship b/w $X$ and $Y$, rather than in subpopulation defined by missing-data pattern.

Selection Model	characterize $y$	missing mechanism
$f(m_i, y_i \vert x_i, \theta, \psi) =$	$f_y(y_i \vert x_i, \theta)$	$\ast f_{m \vert y}(m_i \vert x_i, y_i, \psi)$
$f(m_i, y_i \vert x_i, \xi, \psi) =$	$f(y_i \vert x_i, \xi)$	$\ast f(m_i \vert x_i, \xi)$
Pattern Mixture model	mis 데이터의 다른 패턴들에 의해 정의되는 각각 다른 strata에서의 $y_i \vert x_i$의 분포	probability of different patterns in missingness

missing의 다른 패턴에 따라 $x_i$가 결정이 되고, 그 $x_i$를 기준으로 놓았을 때의 $y_i$의 분포가 궁금.

4.10.5 wk10) Bayesian Model Selection

해당 문제는 prior을 어떻게 고르느냐에 따라서 해결될 수 있음. 이하는 해당 문제에 대한 다양한 해결책들.

4.10.5.1 1. Spike-and-Slab prior

let $X_{n \times p} , Y_{n \times 1}$, then $Y = X \beta$, where $\beta_{p \times 1}$.

p가 많다, 즉 배리어블이 많다는 이야기는 실제 각각의 x의 정보량이 중첩될 가능성이 큼. 그러면 수학적으로는 x’x가 full rank matrix가 아닐 것이며, 이는 곧 몇몇 변수들 간에 서로간의 의존관계가 강하여 의미없는 정보를 포함하는 변수들이 많아질 것. 이러한 의미없는 변수를 삭제하고 실제로 필요한 변수들만 골라내어 y에 대한 inference를 하고 싶음. 이것이 모델 셀렉션 문제이며 이걸 베이지안적으로 풀어내는 것이 곧 Bayesian Model Selection.

무지성 prior로는 $\pi(\beta) \sim N(0, \sigma^2)$가 쓰이지만, 이로는 variable selection이 불가. $\beta$ 중 하나의 component가 0에 가깝게 나왔다고 한들 이것을 0으로 판정할 indicator가 없기 때문. (HPD interval을 구성해서 이것이 0을 포함하면 내치는 식의 방법도 있지만 일단은.)

따라서 다른 prior를 필요로 함. 바로 여기서 사용되는 것이 Spike-and-Slab prior. 이름에서 알 수 있듯이 mixture distribution을 prior로 사용함. $\beta$의 component가 spike 부분에 포함되면 이를 0으로 판정함, 즉 not significant로 판정. 이의 역은 slab part.

이는 곧 prior로 variable selection을 한다는 이야기이다. 즉 이 상황에서는 prior가 패널티로 들어간 것이 된다. 정의적으로 엄밀하게 패널티는 아니지만 사실상의 패널티로 작동. 패널티 term (error penalty)으로 골라내는 것은 full context에서 많이 사용? 이때는 라플라스 프라이어를 쓰고, 노멀을 프라이어 주면? 패널티 텀을 베이지안 인퍼런스로 연결지어서 생각할 수 있지만, 이 배리어블 셀렉션은 디멘션 셀렉션과 연관이 있기 때문에 위와는 정확하게는 다른 개념?

variable selection에는 3가지 방법: 1. 패널티 텀 2. mixture prior 3. 컴퓨테이셔널 (reversible jump, dimension selection) (gradient descent는 아님!)

spike 파트 (아래에서는 $(1-\lambda)N(0, \sigma^2)$) 에는 double exponential을 쓰거나, normal 을 변형해서 사용함. 혹은 극단적으로는 dirac 분포 (point mass) 를 쓸 수도 있음.

이하에서 예시로 제시된 수식은 SS prior이며, 이는 spike와 slab 모두 Normal을 사용하였음.

\[ \pi(\beta) \sim (1-\lambda)N(0, \sigma^2) + \lambda N(0, \omega \sigma^2), \; \; \; \; \; w \gg 1 \]

위에서 $\sigma^2$는 spike variance, $w \sigma^2$는 slab variance.

여기서 $\lambda$가 취할 수 있는 값은 0 아니면 1. 확인할 수 있듯이 1이면 slab part, 즉 significant하고, 0이면 역으로 not significant. 우리는 이에 MCMC 알고리즘을 적용하게 되며, 따라서 MCMC 샘플로 계산을 하면 해당 샘플에서 0인 propotion과 1인 비율이 나오게 될 것. 이때 1인 비율이 0.5 이상이면? 해당 component (변수) 는 significant 하다고 결론짓는 것이 가능하다.

$$ \[\begin{alignat}{3} \pi(\beta, \lambda, \sigma^2, \omega \vert y, x) &\sim \pi(\beta \vert \lambda, \sigma^2, \omega ) \pi(\lambda, \sigma^2, \omega) && f(y \vert x, \beta, \lambda, \sigma^2, \omega) \\ &\sim \pi(\beta \vert \lambda, \sigma^2, \omega ) \pi(\lambda, \sigma^2, \omega) &&L( x, \beta, \lambda, \sigma^2, \omega \vert y) \\ \\ \pi(\lambda) &\sim BETA(1,1), \; \; \; \pi(\sigma^2) \sim \cdot \tag{1}, \; \; \; \omega \sim 1 + GAM(\alpha, \beta) \end{alignat}\] $$

이때 $\sigma^2$는 우리가 임의로 fixed 해서 given으로 잡거나, 위처럼 prior로 해서 시뮬레이션 중에 생산되도록 할수도 있다. 여기선 $\dfrac{1}{U(4,100)}$을 사용.

accept를 하기 위해선 위를 돌리면 됨. 이는 Stochastic Search Variable Selection (SSVS)라고 불림. 이는 GS를 통하여 패러미터를 sequentially update. 이의 결과값은 다음과 같으며, 프로세스는 그 다음과 같다.

\[ (\beta_1 , \cdots, \beta_p, \lambda_1, \cdots, \lambda_p, \sigma^2, \omega) \]

Proceeds:
1. update model parameter $\beta_i^{(t+1)} \sim \pi( \beta_{i} \vert y, x, \beta_{-i}^{(t)}, \lambda_{i}^{(t)}, {\sigma^2}^{(t)}, \omega^{(t)} )$
  - where $\beta_{-i}^{(t)} = \left( \beta_{1}^{(t+1)}, \cdots, \beta_{i-1}^{(t+1)}, \beta_{i+1}^{(t)}, \cdots, \beta_{p}^{(t)} \right)$
  - Simple GS로도 가능하고, MH-within-Gibbs로도 가능함
2. update $\lambda_I^{(t+1)} \sim \pi(\lambda_i \vert y, x, \lambda_{-i}^{(t)}, \beta_{i}^{(t+1)}, {\sigma^2}^{(t)}, \omega^{(t)} )$
  - where
    $P( {i}^{(t+1)} = 1 y, x, {-i}^{(t)}, _{i}^{(t+1)}, {^2}{(t)}, ^{(t)} )= {a+b} BER( {a+b} ) $.
    - $a = \pi( \beta_{i}^{(t+1)} \vert y, x, \lambda_{i}^{(t+1)}=1, {\sigma^2}^{(t)}, \omega^{(t)} )$.
    - $b = \pi( \beta_{i}^{(t+1)} \vert y, x, \lambda_{i}^{(t+1)}=0, {\sigma^2}^{(t)}, \omega^{(t)} )$,
3. update $\sigma^2$
4. update $\omega$

4.10.5.2 2. Horseshoe prior (Scale-mixture prior)

distribution에서 scale이란 Variance.

위와 동일 케이스 가정. 그 경우

$$ \[\begin{alignat}{4} \pi(\beta \vert y, x) \propto f(y \vert x, \beta) \pi(\beta), \; \; \; \; \; &\pi(\beta) && \sim N(0, \sigma^2) \\ \Longrightarrow &\pi(\beta_i \vert \tau, \lambda_i) && \sim N(0, \tau^2 \lambda_i^2) \end{alignat}\] $$

where $pi(\tau^2), \pi(\lambda_i^2) \sim Cauchy^{+}(0,1)$

이때 common variance component $\tau$는 각 component마다 공유하는 1개의 variance component, 그리고 각 component마다 indiviually 고유한 individual parameter variance component $\lambda_i$를 설정한 것.

4.10.6 Autologistic model

4.10.7 wk10) Bayesian Model Averaging

해당 상황에서 연구자는 다양한 모델 예측 후보를 생각해볼 수 있음. 보통은 프로세스를 거쳐 이 모델들 중의 하나를 선택하게 됨. 하지만 완벽한 모델이라는 건 (보통) 존재할 수 없음, 어떤 모델 후보를 선택하든 해당 후보가 내포하고 있는 uncertainty가 존재하며 이를 수용하게 됨. 따라서 모델을 선택한다는 것은 동시에 over-confidence inference 문제를 발생시킨다. 따라서 모델 후보군을 하나만 골라야 한다는 고정관념을 벗어나 다양한 모델 후보군들 각각을 동시에 반영하자. 이 동시 반영할 때 각 모델이 내포하고 있는 확률 (uncertainty)에 의해 각 모델의 반영 정도를 가감하게 된다.

BMA는 패러미터 estimate를 획득할 때, 이러한 model uncertainty를 설명하기 위한 일관된 메커니즘을 제공한다.

given 데이터 $D$, posterior prob of $\mathcal{M}_k$ $= P(\mathcal{M}_k \vert D) = \dfrac{L(D \vert \mathcal{M}_k) P(\mathcal{M}_k)}{\sum_{k=1}^k L(D \vert \mathcal{M}_k) P(\mathcal{M}_k)}$.

이때 marginal likelihood under $\mathcal{M}_k)$ $L(D \vert \mathcal{M}_k) = \int L(D \vert \theta \mathcal{M}_k) \pi(\theta \vert \mathcal{M}_k) d \theta$이며, integral 안의 수식은 posterior of model의 상수배

In brief, BMA는 model uncertainty를 설명할 수 있는 posterior density를 획득하기 위해 integral을 취한다 (model에 대해 마지널化). (각각의 모델에 대한 model uncertainty를 구한다)

이를 통해 최종적으로 posterior sample 같은 경우에는 각 모델 별로 model probability에 posterior sample의 probability를 다 더해준 값이 실제로 우리의 $\theta$에 대한 post가 된다.

즉슨 BMA란 다양한 모델 후보군들이 존재할 때, 그 어떤 상황에서도 robust inference를 가능케 하는 tool이 바로 BMA.

4.10.7.1 Ex: BMA-CRM

CRM 시작 전에 skeleton 정하고 시작하는 건 자명. 근데 이 skeleton이 잘못 선정되었다면 제대로된 도즈 selection이 불가능해지므로 skeleton의 선정이 잘못되어 있다면 이는 치명적임. 상식적으로, 하나의 skeleton으로만 도즈 셀렉션을 진행한다면 문제가 생길 확률은 당연히 높음. 이런 리스크를 희석하기 위해 skeleton을 다수를 정하고 CRM을 시작하면 이런 문제를 다소 회피할 수 있지 않을까? 이때 이 각각의 스켈레톤 하나하나를 모델로 인식한다. 이 각각의 모델에 따라서 CRM을, 즉 도즈레벨을 업데이트할 확률을, 즉 업데이트 할 때 패러미터 evaluation을 하는데, 그때 나오는 패러미터 값과 그 각각의 모델 probability를 비교하여 그 모델 averaging을 해주면 그 어떤 상황이 와도 굉장히 robust 한 값을 획득할 수 있을 것.

Main research question
Justification for your research question (why is it important to answer the question?)
Data source
Data analysis
Summary of the data analysis results and conclusion
Appendix (if needed) – R scripts (scripts or codes for any other software)

1 – Technical details regarding the statistical tools used in the analysis

비교적 오랜 기간 데이터가 잘 정립된 MLB 기록 활용을 위해 수업 시 활용하였던 Lahman package(R) 사용

야구는 공격과 수비로 이루어지며, 따라서 분석을 진행함에 있어

야구는

야구 스탯들 간에 상당한 수준의 선형성이 보장되어 이

타자의 가치 = α∗Batting+β∗Fielding, α,β는 임의의 패러미터

선수의 타격능력은 이른바 클래식 스탯으로 불리는 다양한 구형 통계량으로도 표기하는데 문제가 없지만, 수비능력은 선수별로 할당된 수비범위가 천차만별이며 선수가 수비시도를 하지 않으면 선수의 실책으로 이어지지 않는다는 점 때문에 선수 개별의 수비능력이 객관적으로 평가되기 시작한 것은 구장의 정보를 훨씬 자세하게 담을 수 있게 된 2000년대 중반부 이후부터의 이야기.

최신야구에서 선수 수비능력의 평가는 주로 Ultimate Zone Rating (UZR) 로 이루어지며, 해당 스탯은 ARM (달린거리), DPR (병살), RngR (수비커버리지), ErrR (에러빈도) 등 수비에 관련된 스탯을 총집합하여 망라하는 고밀도 스탯이다. 그러나 해당 스탯의 계산은 2002년 BIS (Baseball Info Solutions)라는 회사에서 제공하는 유료 데이터셋과 15년 도입된 스탯캐스트 데이터에 거의 전적으로 의존하고 있다. 스탯캐스트 데이터는 민간에 어느정도는 공개되어 있어 접근이 불가능하지 않지만 (https://baseballsavant.mlb.com/statcast_leaderboard), BIS 데이터는 접근이 어렵다.

이와 같은 이유로 선수별 수비 스탯을 구하기를 시도하기보단 팀별 수비력에 대한 척도인 Defensive Efficiency Ratio (DER)를 사용하고자 한다. 최대 12개의 팀인만큼 팀 간의 차이를 포착하기 쉬우며, 12개의 팀으로 표준화되니만큼 아웃라이어들이 평준화되어 전반적인 경향성으로 기능하는 것을 기대해볼 수 있다. DER의 수식은 다음과 같다.

DER의 계산법은 이하와 같다.

DER=1−(Hits+Reached.On.Error−HomeRunsPlate.Appearance−BB(Walks)−Strike.Out−Hit.By.Pitch−HomeRuns)

Teams %>% mutate(., DER = 1-((H + E - HR)/((AB + SF) - SO - HR))) %>% ##select(., yearID, teamID, Rank, SO, SOA) select(.,yearID, teamID, franchID, Rank, G,DER) -> Teams_DER

물론 같은 팀에 속했다는 이유만으로 모든 선수들에게 동급의 수비스탯을 배정하는 건 합리적이라고 말하기 어렵다. 팀의 수비에 기여하는 정도가 높은 선수가 있다면 낮은 선수도 있을 것이 자명하기 때문이다. 따라서 팀별로 획득한 DER을 수비에 대한 클래식 스탯인 각 선수의 Fielding Percentage(FPCT) 나 Range Factor(RF)의 비율로 스케일링해서 부여하자. 두 스탯은 각각 수비능력과 개인의 수비범위 평가를 위해 시도되었던 스탯들이지만, 전자는 개인의 수비범위가 좁으면 더 좋은 값이 나온다는 한계, 후자는 공이 본인 위치로 떨어졌을 때 스탯계산에서 이득을 본다는 한계를 넘지 못해 좋은 스탯으로는 평가받지 못했던 값들이다. 그러나 팀 단위로 한번 수비력을 표준화한 후 팀 내에서 상대적인 기여도를 보는 식으로 보정이 한 번 들어갔으므로 어느정도 기준선으로서는 기능할 것이라고 기대된다.

이를 위해 선수생활 중의 메인 수비포지션 지정하고 해당 포지션에서의 통계량만 사용.

Selection Model	characterize \(y\)	missing mechanism
\(f(m_i, y_i \vert x_i, \theta, \psi) =\)	\(f_y(y_i \vert x_i, \theta)\)	\(\ast f_{m \vert y}(m_i \vert x_i, y_i, \psi)\)
\(f(m_i, y_i \vert x_i, \xi, \psi) =\)	\(f(y_i \vert x_i, \xi)\)	\(\ast f(m_i \vert x_i, \xi)\)
Pattern Mixture model	mis 데이터의 다른 패턴들에 의해 정의되는 각각 다른 strata에서의 \(y_i \vert x_i\)의 분포	probability of different patterns in missingness