s trts r - PowerPoint PPT Presentation

❙❦❡✇✲◆♦r♠❛❧ ❉❛t❛ ❙✉♣♣♦s❡ ✇❡ ❣❡♥❡r❛t❡ ❞❛t❛ ❢r♦♠ ❛ s❦❡✇✲♥♦r♠❛❧ ❞✐str✐❜✉t✐♦♥ ✶ ❙◆ ( µ, ω ✷ , α ) ✱ ✇❤❡r❡ µ ∈ ℜ ✱ ω > ✵✱ ❛♥❞ α ∈ ℜ ❛r❡ ❧♦❝❛t✐♦♥✱ s❝❛❧❡✱ ❛♥❞ s❦❡✇♥❡ss ♣❛r❛♠❡t❡rs α > ✵ ⇒ ♣♦s✐t✐✈❡ s❦❡✇♥❡ss✱ α < ✵ ⇒ ♥❡❣❛t✐✈❡ s❦❡✇♥❡ss✱ ❛♥❞ ✇❤❡♥ α = ✵✱ t❤❡ ❞❡♥s✐t② r❡❞✉❝❡s t♦ ❛ s②♠♠❡tr✐❝ ♥♦r♠❛❧ ❞✐str✐❜✉t✐♦♥ ❋♦r ❞❡t❛✐❧s ♦♥ ❇❛②❡s✐❛♥ ❛♣♣r♦❛❝❤❡s t♦ ✜tt✐♥❣ ❙◆ ♠♦❞❡❧s✱ s❡❡ ❋rü❤✇✐rt❤✲❙❝❤♥❛tt❡r ❛♥❞ P②♥❡ ✭✷✵✶✵✮ ❛♥❞ ◆❡❡❧♦♥ ❡t ❛❧✳ ✭✷✵✶✺✮ ❋♦r ♥♦✇✱ s✉♣♣♦s❡ ✇❡ ✐❣♥♦r❡ s❦❡✇♥❡ss ❛♥❞ ✜t ❛♥ ♦r❞✐♥❛r② ❧✐♥❡❛r r❡❣r❡ss✐♦♥ t♦ t❤❡ ❞❛t❛ ❙❡❡ ❘❡s✐❞✉❛❧ ❉✐❛❣♥♦st✐❝s ✇✐t❤ ❙◆ ❉❛t❛✳r ❢♦r ❞❡t❛✐❧s ✶ ❖✬❍❛❣❛♥ ❛♥❞ ▲❡♦♥❛r❞ ✭✶✾✼✻✮❀ ❆③③❛❧✐♥✐✱ ✶✾✽✺ ✶✶ ✴ ✶✻✵

P❧♦ts ♦❢ ❙t❛♥❞❛r❞✐③❡❞ ❘❡s✐❞✉❛❧s True Errors (e) Density of Standardized Residuals α = −5 MLE 0.4 0.15 ω = 4 Bayes 0.3 0.10 Density Density 0.2 0.05 0.1 0.00 0.0 −15 −10 −5 0 −4 −2 0 2 y Standardized Residual QQ Plot of Standardized Residuals 3 MLE ● ● Bayes 2 ● Observed Quantile ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 ● ● ●● ●● −3 ● ● −3 −2 −1 0 1 2 3 Normal Quantile ✶✷ ✴ ✶✻✵

Pr♦❜✐t ❛♥❞ ▲♦❣✐t ▼♦❞❡❧s ❢♦r ❇✐♥❛r② ❉❛t❛ ∗ ❈♦♥s✐❞❡r t❤❡ ❢♦❧❧♦✇✐♥❣ ♣r♦❜✐t ♠♦❞❡❧ ❢♦r ❛ ❞✐❝❤♦t♦♠♦✉s ♦✉t❝♦♠❡ Y i ✿ Φ − ✶ [ Pr ( Y i = ✶ )] ① T = i β , i = ✶ , . . . , n , ✇❤❡r❡ Φ( · ) ❞❡♥♦t❡s t❤❡ ❈❉❋ ♦❢ ❛ st❛♥❞❛r❞ ♥♦r♠❛❧ r❛♥❞♦♠ ✈❛r✐❛❜❧❡ ❲❡ ❝❛♥ r❡♣r❡s❡♥t t❤❡ ♠♦❞❡❧ ✈✐s✲à✲✈✐s ❛ ❧❛t❡♥t ✈❛r✐❛❜❧❡ Z i s✉❝❤ t❤❛t ◆ ( ① T ∼ i β , ✶ ) ❛♥❞ Z i Y i = ✶ ⇐ ⇒ Z i > ✵ ✐♠♣❧②✐♥❣ t❤❛t Pr ( Z i > ✵ ) = Φ( ① T Pr ( Y i = ✶ ) = i β ) ✶✸ ✴ ✶✻✵

▲❛t❡♥t ❱❛r✐❛❜❧❡ ■♥t❡r♣r❡t❛t✐♦♥ ♦❢ Pr♦❜✐t ▼♦❞❡❧ 0.4 0.3 f(z i ) 0.2 Pr(Z i > 0) = Pr(Y i = 1) 0.1 0.0 T β x i −2 0 2 4 z i ✶✹ ✴ ✶✻✵

❆❧❜❡rt ❛♥❞ ❈❤✐❜ ✭✶✾✾✸✮ ❉❛t❛✲❆✉❣♠❡♥t❡❞ ❙❛♠♣❧❡r ❆❧❜❡rt ❛♥❞ ❈❤✐❜ ✭✶✾✾✸✮ t❛❦❡ ❛❞✈❛♥t❛❣❡ ♦❢ t❤✐s ❧❛t❡♥t ✈❛r✐❛❜❧❡ str✉❝t✉r❡ t♦ ❞❡✈❡❧♦♣ ❛♥ ❡✣❝✐❡♥t ❞❛t❛✲❛✉❣♠❡♥t❡❞ ●✐❜❜s s❛♠♣❧❡r ❢♦r ♣r♦❜✐t r❡❣r❡ss✐♦♥ ❉❛t❛ ❛✉❣♠❡♥t❛t✐♦♥ ✐s ❛ ♠❡t❤♦❞ ❜② ✇❤✐❝❤ ✇❡ ✐♥tr♦❞✉❝❡ ❛❞❞✐t✐♦♥❛❧ ✭♦r ✏❛✉❣♠❡♥t❡❞✑✮ ✈❛r✐❛❜❧❡s✱ ❩ = ( Z ✶ , . . . , Z n ) T ✱ ❛s ♣❛rt ♦❢ t❤❡ ●✐❜❜s s❛♠♣❧❡r t♦ ❢❛❝✐❧✐t❛t❡ s❛♠♣❧✐♥❣ ❉❛t❛ ❛✉❣♠❡♥t❛t✐♦♥ ✐s ✉s❡❢✉❧ ✇❤❡♥ t❤❡ ❝♦♥❞✐t✐♦♥❛❧ ❞❡♥s✐t② π ( β | ② ) ✐s ✐♥tr❛❝t❛❜❧❡✱ ❜✉t t❤❡ ❥♦✐♥t ♣♦st❡r✐♦r π ( β , ③ | ② ) ✐s ❡❛s② t♦ s❛♠♣❧❡ ❢r♦♠ ✈✐❛ ●✐❜❜s✱ ✇❤❡r❡ ③ ✐s ❛♥ n × ✶ ✈❡❝t♦r ♦❢ r❡❛❧✐③❛t✐♦♥s ♦❢ ❩ ✶✺ ✴ ✶✻✵

❉❛t❛ ❆✉❣♠❡♥t❛t✐♦♥ ❙❛♠♣❧❡r ■♥ ♣❛rt✐❝✉❧❛r✱ s✉♣♣♦s❡ ✐t✬s str❛✐❣❤t❢♦r✇❛r❞ t♦ s❛♠♣❧❡ ❢r♦♠ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧s π ( β | ② , ③ ) ❛♥❞ π ( ③ | ② , β ) ❚❤❡♥ ✇❡ ❝❛♥ ❛♣♣❧② ●✐❜❜s s❛♠♣❧✐♥❣ t♦ ♦❜t❛✐♥ t❤❡ ❥♦✐♥t ♣♦st❡r✐♦r π ( β , ③ | ② ) ❆❢t❡r ❝♦♥✈❡r❣❡♥❝❡✱ t❤❡ s❛♠♣❧❡s ♦❢ β ✱ { β ( ✶ ) , . . . , β ( T ) } ✱ ✇✐❧❧ ❝♦♥st✐t✉t❡ ❛ ▼♦♥t❡ ❈❛r❧♦ s❛♠♣❧❡ ❢r♦♠ π ( β | ② ) ◆♦t❡ t❤❛t ✐❢ β ❛♥❞ ❨ ❛r❡ ❝♦♥❞✐t✐♦♥❛❧❧② ✐♥❞❡♣❡♥❞❡♥t ❣✐✈❡♥ ❩ ✱ s♦ t❤❛t π ( β | ② , ③ ) = π ( β | ③ ) ✱ t❤❡♥ t❤❡ s❛♠♣❧❡r ♣r♦❝❡❡❞s ✐♥ t✇♦ st❛❣❡s✿ ✶ ❉r❛✇ ③ ❢r♦♠ π ( ③ | ② , β ) ✷ ❉r❛✇ β ❢r♦♠ π ( β | ③ ) ✶✻ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r Pr♦❜✐t ▼♦❞❡❧ ∗ ❚❤❡ ❞❛t❛ ❛✉❣♠❡♥t❡❞ s❛♠♣❧❡r ♣r♦♣♦s❡❞ ❜② ❆❧❜❡rt ❛♥❞ ❈❤✐❜ ♣r♦❝❡❡❞s ❜② � � β ✵ , ❚ − ✶ ❛ss✐❣♥✐♥❣ ❛ ◆ p ♣r✐♦r t♦ β ❛♥❞ ❞❡✜♥✐♥❣ t❤❡ ♣♦st❡r✐♦r ✈❛r✐❛♥❝❡ ✵ � � − ✶ ❚ ✵ + ❳ T ❳ ♦❢ β ❛s ❱ = ◆♦t❡ t❤❛t ❜❡❝❛✉s❡ ❱❛r ( Z i ) = ✶✱ ✇❡ ❝❛♥ ❞❡✜♥❡ ❱ ♦✉ts✐❞❡ t❤❡ ●✐❜❜s ❧♦♦♣ ◆❡①t✱ ✇❡ ✐t❡r❛t❡ t❤r♦✉❣❤ t❤❡ ❢♦❧❧♦✇✐♥❣ ●✐❜❜s st❡♣s✿ ✶ ❋♦r i = ✶ , . . . , n ✱ s❛♠♣❧❡ z i ❢r♦♠ ❛ ◆ ( ① T i β , ✶ ) ❞✐str✐❜✉t✐♦♥ tr✉♥❝❛t❡❞ ❜❡❧♦✇ ✭❛❜♦✈❡✮ ❜② ✵ ❢♦r y i = ✶ ( y i = ✵ ) � � ❚ ✵ β ✵ + ❳ T ③ ✷ ❙❛♠♣❧❡ β ❢r♦♠ ◆ p ( ♠ , ❱ ) ✱ ✇❤❡r❡ ♠ = ❱ ❛♥❞ ❱ ✐s ❞❡✜♥❡❞ ❛❜♦✈❡ ◆♦t❡✿ ❈♦♥❞✐t✐♦♥❛❧ ♦♥ ❩ ✱ β ✐s ✐♥❞❡♣❡♥❞❡♥t ♦❢ ❨ ✱ s♦ ✇❡ ❝❛♥ ✇♦r❦ s♦❧❡❧② ✇✐t❤ t❤❡ ❛✉❣♠❡♥t❡❞ ❧✐❦❡❧✐❤♦♦❞ ✇❤❡♥ ✉♣❞❛t✐♥❣ β ❙❡❡ Pr♦❜✐t✳r ❢♦r ❞❡t❛✐❧s ✶✼ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r Pr♦❜✐t ●✐❜❜s ❙❛♠♣❧❡r Gibbs for Probit Regression Model # Priors beta0<- rep (0,p) # Prior mean of beta (of dimension p) T0<- diag (.01,p) # Prior precision of beta # Inits beta<- rep (0,p) z<- rep (0,n) # Latent normal variables ns<- table (y) # Category sample sizes # MCMC info analogous to linear reg. code # Posterior var of beta -- Note: can calculate outside of loop vbeta<- solve (T0+ crossprod (X,X)) ######### # Gibbs # ######### tmp<- proc.time () # Store current time for (i in 1:nsim){ # Update latent normals, z, from truncated normal using inverse-CDF method muz<-X%*%beta z[y==0]<- qnorm ( runif (ns[1],0, pnorm (0,muz[y==0],1)),muz[y==0],1) z[y==1]<- qnorm ( runif (ns[2], pnorm (0,muz[y==1],1),1),muz[y==1],1) # Alternatively, can use rtnorm function from msm package -- this is slower # z[y==0]<-rtnorm(n0,muz[y==0],1,-Inf,0) # z[y==1]<-rtnorm(n1,muz[y==1],1,0,Inf) # Update beta mbeta <- vbeta%*%(T0%*%beta0+ crossprod (X,z)) # Posterior mean of beta beta<- c ( rmvnorm (1,mbeta,vbeta)) ################# # Store Results # ################# if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta } if (i%%100==0) print (i) } proc.time ()-tmp # MCMC run time -- 0.64 seconds to run 1000 iterations with n=1000 ✶✽ ✴ ✶✻✵

❊①❛♠♣❧❡ ❚❤❡ ♣r♦❣r❛♠ Pr♦❜✐t✳r ✜ts t❤❡ ❢♦❧❧♦✇✐♥❣ ♣r♦❜✐t ♠♦❞❡❧✿ ❇❡r♥ ( π i ) Y i ∼ Φ − ✶ ( π i ) = β ✵ + β ✶ x i , i = ✶ , . . . , ✶✵✵✵ . ❚❤❡ r❡s✉❧ts✱ ❜❛s❡❞ ♦♥ ✶✵✵✵ ✐t❡r❛t✐♦♥s ✇✐t❤ ❛ ❜✉r♥✲✐♥ ♦❢ ✺✵✵✱ ❛r❡✿ ❚❛❜❧❡ ✷✿ ❘❡s✉❧ts ❢♦r Pr♦❜✐t ▼♦❞❡❧✳ ❚r✉❡ P♦st❡r✐♦r ▼❡❛♥ ✭❙❉✮ † P❛r❛♠❡t❡r ❱❛❧✉❡ ▼▲❊ ✭❙❊✮ β ✵ − ✵ . ✺ − ✵ . ✻✹ ( ✵ . ✵✼ ) − ✵ . ✻✹ ( ✵ . ✵✼ ) β ✶ ✵ . ✺ ✵ . ✺✺ ( ✵ . ✵✹ ) ✵ . ✺✻ ( ✵ . ✵✹ ) † ❇❛s❡❞ ♦♥ ❆❧❜❡rt ❛♥❞ ❈❤✐❜ ❞❛t❛ ❛✉❣♠❡♥t❛t✐♦♥ s❛♠♣❧❡r✳ ✶✾ ✴ ✶✻✵

t ✲▲✐♥❦ ▼♦❞❡❧s ❆❧❜❡rt ❛♥❞ ❈❤✐❜ ❛❧s♦ ❞✐s❝✉ss ❡①t❡♥s✐♦♥s t♦ s♦✲❝❛❧❧❡❞ t ✲❧✐♥❦ ♠♦❞❡❧s t❤❛t ❛❧❧♦✇ Z t♦ ❛ss✉♠❡ ❛ t ✲❞✐str✐❜✉t✐♦♥ ✇✐t❤ ❤❡❛✈✐❡r t❛✐❧s t❤❛♥ t❤❡ ♥♦r♠❛❧ ❞✐str✐❜✉t✐♦♥ ❚❤❡ ❝❤♦✐❝❡ ♦❢ ❞❢ ✇✐❧❧ ❝♦rr❡s♣♦♥❞ t♦ ❞✐✛❡r❡♥t ❧✐♥❦ ❢✉♥❝t✐♦♥s ❆ t ✶ ❞✐str✐❜✉t✐♦♥ ❝♦rr❡s♣♦♥❞s t♦ ❛ ❈❛✉❝❤② ❧✐♥❦ ❆ t ✽ ❞✐str✐❜✉t✐♦♥ ❛♣♣r♦①✐♠❛t❡s ❛ ✭s❝❛❧❡❞✮ ❧♦❣✐st✐❝ ❞✐str✐❜✉t✐♦♥ ❆♥❞✱ ✐♥ t❤❡ ❧✐♠✐t✱ t ∞ ✐♠♣❧✐❡s ❛ ♣r♦❜✐t ❧✐♥❦ ❇② ✈❛r②✐♥❣ t❤❡ ❞❢s✱ ✇❡ ♣❡r♠✐t ✢❡①✐❜✐❧✐t② ✐♥ t❤❡ ❝❤♦✐❝❡ ♦❢ ❧✐♥❦ ❢✉♥❝t✐♦♥ ✷✵ ✴ ✶✻✵

▲❛t❡♥t ❱❛r✐❛❜❧❡ ■♥t❡r♣r❡t❛t✐♦♥ ♦❢ t ✲▲✐♥❦ ▼♦❞❡❧ 0.4 T β ) t 3 (x i T β ,1) N(x i 0.3 f(z i ) 0.2 0.1 Pr(Z i > 0) = Pr(Y i = 1) 0.0 0 T β x i z i ✷✶ ✴ ✶✻✵

t ✲▲✐♥❦ ✈s ▲♦❣✐st✐❝ ◗✉❛♥t✐❧❡s ■♥ ♣❛rt✐❝✉❧❛r✱ ❛ t ✽ / ✵ . ✻✸✹ ❛♣♣r♦①✐♠❛t❡s ❛ ❧♦❣✐st✐❝ ❞✐str✐❜✉t✐♦♥ t 3 t 8 t 16 Logistic*0.634 5 logistic quantile 0 −5 −4 −2 0 2 4 t quantile ✷✷ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r t❤❡ t ✲▲✐♥❦ ▼♦❞❡❧ ∗ ❘❡❝❛❧❧ t❤❛t ❛ t ✲❞✐str✐❜✉t✐♦♥ ❛r✐s❡s ❛s ❛ s❝❛❧❡❞ ♠✐①t✉r❡ ♦❢ ♥♦r♠❛❧s✱ ✇❤❡r❡ t❤❡ s❝❛❧❡ ✭✐✳❡✳✱ ✈❛r✐❛♥❝❡✮ ♣❛r❛♠❡t❡r ❢♦❧❧♦✇s ❛♥ ■● ❞✐str✐❜✉t✐♦♥ ❊q✉✐✈❛❧❡♥t❧②✱ t❤❡ ♣r❡❝✐s✐♦♥ ♦❢ t❤❡ ♥♦r♠❛❧ ✈❛r✐❛❜❧❡s ✐s ❣❛♠♠❛ ■♥ ♣❛rt✐❝✉❧❛r✱ ❛ t ν r❛♥❞♦♠ ✈❛r✐❛❜❧❡✱ Z i ( i = ✶ , . . . , n ) ✱ ❝❛♥ ❜❡ ❣❡♥❡r❛t❡❞ ❜② ✶ ●❡♥❡r❛t✐♥❣ λ i ❢r♦♠ ●❛ ( ν/ ✷ , ν/ ✷ ) ✭♣❛r❛♠❡t❡r✐③❡❞ ❛s r❛t❡ ✮ � � i β , λ − ✶ ① T ✷ ●❡♥❡r❛t✐♥❣ Z i | λ i ❢r♦♠ ◆ i ▼❛r❣✐♥❛❧❧②✱ Z i ✐s ❞✐str✐❜✉t❡❞ ❛s t ν ( ① T i β ) ❛♥❞ t❤❡ ♦r✐❣✐♥❛❧ ❜✐♥❛r② Y i ✐s ♠♦❞❡❧❡❞ ✉s✐♥❣ t❤❡ ❝♦rr❡s♣♦♥❞✐♥❣ t ν ✲❧✐♥❦ ❢✉♥❝t✐♦♥ ✷✸ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r t❤❡ t ✲▲✐♥❦ ▼♦❞❡❧ ∗ ❚❤✐s ❧❡❛❞s t♦ t❤❡ ❢♦❧❧♦✇✐♥❣ ♠♦❞✐✜❝❛t✐♦♥ ❢♦r t❤❡ ♣r♦❜✐t ●✐❜❜s s❛♠♣❧❡r ✶ ❋♦r i = ✶ , . . . , n ✱ s❛♠♣❧❡ z i ❢r♦♠ ✐ts tr✉♥❝❛t❡❞ ◆ ( ① T i β , λ − ✶ i ) ❞✐str✐❜✉t✐♦♥ ❛♥❛❧♦❣♦✉s t♦ t❤❡ ❡❛r❧✐❡r s❛♠♣❧❡r ✷ ❙❛♠♣❧❡ β ❢r♦♠ ◆ p ( ♠ , ❱ ) ✱ ✇❤❡r❡ � � − ✶ ❚ ✵ + ❳ T ❲ ❳ ❱ = � � ❚ ✵ β ✵ + ❳ T ❲ ③ ♠ = ❱ ✭❛♥❛❧♦❣♦✉s t♦ ❲▲❙✮ ❲ = ❞✐❛❣ ( λ i ) � � ( ν + ✶ ) / ✷ , ( ν + ( z i − ① T i β ) ✷ ) / ✷ ✸ ❋♦r i = ✶ , . . . , n ✱ s❛♠♣❧❡ λ i ❢r♦♠ ●❛ ✹ ❖♣t✐♦♥❛❧❧②✱ ♣❧❛❝❡ ❛ ❞✐s❝r❡t❡ ✉♥✐❢♦r♠ ♣r✐♦r ♦♥ ν ❛♥❞ ❛♣♣❧② t❤❡ ❞✐s❝r❡t❡ ✈❡rs✐♦♥ ♦❢ ❇❛②❡s✬ ❚❤❡♦r❡♠ t♦ ❣❡t t❤❡ ♣♦st❡r✐♦r ♣r♦❜❛❜✐❧✐t✐❡s ✺ ❋♦r t❤❡ ❧♦❣✐st✐❝ ❛♣♣r♦①✐♠❛t✐♦♥✱ s❡t ν = ✽ ❛♥❞ ❞✐✈✐❞❡ t❤❡ ♣♦st❡r✐♦r ♠❡❛♥s ❛♥❞ ❙❉s ❜② ✵✳✻✸✹ t♦ r❡❝♦✈❡r t❤❡ ❧♦❣✐st✐❝ r❡s✉❧ts ❙❡❡ t ✲▲✐♥❦ ▲♦❣✐t✳r ❢♦r ❞❡t❛✐❧s ✷✹ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r t ✲▲✐♥❦ ▲♦❣✐st✐❝ ▼♦❞❡❧ Gibbs Sampler for t-Link Logit Model # A&C assign diffuse (improper) priors for beta, so no explicit specification # Initial Values lambda<- rep (1,n) # Weights nu<-8 # t df (8 ~ logistic) -- assume fixed z<- rep (0,n) # Latent z vector ns<- table (y) # Category sample szies ########### # Gibbs # ########### tmp<- proc.time () # Store current time for (i in 1:nsim) { # Update z using inverse-CDF method muz<-X%*%beta z[y==0]<- qnorm ( runif (ns[1],0, pnorm (0,muz[y==0], sqrt (1/lambda[y==0]))),muz[y==0], sqrt (1/lambda[y==0])) z[y==1]<- qnorm ( runif (ns[2], pnorm (0,muz[y==1], sqrt (1/lambda[y==1])),1),muz[y==1], sqrt (1/lambda[y==1])) vbeta<- solve ( crossprod ( sqrt (lambda)*X)) # Can no longer update outside Gibbs loop betahat<-vbeta%*%( crossprod (lambda*X,z)) beta<- c ( rmvnorm (1,betahat,vbeta)) # Update lambda lambda<- rgamma (n,(nu+1)/2,(nu+(z-X%*%beta)^2)/2) ################# # Store Results # ################# if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta } if (i%%100==0) print (i) } proc.time ()-tmp # MCMC run time -- 1 sec to run 1000 iterations with n=1000 # Results mbeta<- colMeans (Beta/.634) # Correction factor is 1/.634 sbeta<- apply (Beta/.634,2,sd) ✷✺ ✴ ✶✻✵

❊①❛♠♣❧❡ ❚❤❡ ♣r♦❣r❛♠ t ✲▲✐♥❦ ▲♦❣✐t✳r ❣❡♥❡r❛t❡s ✶✵✵✵ ♦❜s❡r✈❛t✐♦♥s ❢r♦♠ t❤❡ ❢♦❧❧♦✇✐♥❣ ❧♦❣✐st✐❝ ♠♦❞❡❧✿ ∼ ❇❡r♥ ( π i ) Y i ❧♦❣✐t ( π i ) = β ✵ + β ✶ x i , i = ✶ , . . . , ✶✵✵✵ . ❚❤❡ r❡s✉❧ts✱ ❜❛s❡❞ ♦♥ ✶✵✵✵ ✐t❡r❛t✐♦♥s ✇✐t❤ ❛ ❜✉r♥✲✐♥ ♦❢ ✺✵✵✱ ❛r❡✿ ❚❛❜❧❡ ✸✿ ❘❡s✉❧ts ❢♦r ▲♦❣✐st✐❝ ▼♦❞❡❧ ✇✐t❤ t ✲▲✐♥❦✳ ❚r✉❡ P♦st❡r✐♦r ▼❡❛♥ ✭❙❉✮ † P❛r❛♠❡t❡r ❱❛❧✉❡ ▼▲❊ ✭❙❊✮ β ✵ − ✶ − ✶ . ✵✽ ( ✵ . ✵✽ ) − ✶ . ✵✽ ( ✵ . ✵✽ ) β ✶ ✶ ✵ . ✾✷ ( ✵ . ✵✾ ) ✵ . ✾✸ ( ✵ . ✵✾ ) † ❇❛s❡❞ ♦♥ ❆❧❜❡rt ❛♥❞ ❈❤✐❜ t ✲❧✐♥❦ ♠♦❞❡❧✳ ✷✻ ✴ ✶✻✵

▲♦❣✐st✐❝ ❘❡❣r❡ss✐♦♥ ❯s✐♥❣ Pó❧②❛✲●❛♠♠❛ ▲❛t❡♥t ❱❛r✐❛❜❧❡s P♦❧s♦♥ ❡t ❛❧✳ ✭✷✵✶✷✮ ♣r♦♣♦s❡❞ ❛ ❛❧t❡r♥❛t✐✈❡ ●✐❜❜s s❛♠♣❧❡r ❢♦r ❧♦❣✐st✐❝ ❛♥❞ ♥❡❣❛t✐✈❡ ❜✐♥♦♠✐❛❧ ♠♦❞❡❧s ❚❤❡ ❛♣♣r♦❛❝❤ ✐♥tr♦❞✉❝❡s ❛ ✈❡❝t♦r ♦❢ ❧❛t❡♥t ✈❛r✐❛❜❧❡s✱ Z i ✱ t❤❛t ❛r❡ s❝❛❧❡ ♠✐①t✉r❡s ♦❢ ♥♦r♠❛❧s ✇✐t❤ ✐♥❞❡♣❡♥❞❡♥t Pó❧②❛✲●❛♠♠❛ ♣r❡❝✐s✐♦♥ t❡r♠s r❛t❤❡r t❤❛♥ ●❛♠♠❛ ♣r❡❝✐s✐♦♥ t❡r♠s ❛s ✐♥ t❤❡ t ✲❧✐♥❦ ♠♦❞❡❧ ❆ r❛♥❞♦♠ ✈❛r✐❛❜❧❡ ω ✐s s❛✐❞ t♦ ❤❛✈❡ ❛ Pó❧②❛✲●❛♠♠❛ ❞✐str✐❜✉t✐♦♥ ✇✐t❤ ♣❛r❛♠❡t❡rs b > ✵ ❛♥❞ c ∈ ℜ ✱ ✐❢ � ∞ ✶ g k d ω ∼ P● ( b , c ) = ( k − ✶ / ✷ ) ✷ + c ✷ / ( ✹ π ✷ ) , ✷ π ✷ k = ✶ ✇❤❡r❡ t❤❡ g k ✬s ❛r❡ ✐♥❞❡♣❡♥❞❡♥t❧② ❞✐str✐❜✉t❡❞ ❛❝❝♦r❞✐♥❣ t♦ ❛ ●❛ ( b , ✶ ) ❞✐str✐❜✉t✐♦♥ ✷✼ ✴ ✶✻✵

Pó❧②❛✲●❛♠♠❛ ❉❡♥s✐t② P❧♦t 4 PG(1, 0) Ga(2, 10) 3 f(x) 2 1 0 0.0 0.5 1.0 1.5 x ✷✽ ✴ ✶✻✵

Pr♦♣❡rt✐❡s ♦❢ t❤❡ Pó❧②❛✲●❛♠♠❛ ❉✐str✐❜✉t✐♦♥ P♦❧s♦♥ ❡t ❛❧✳ ❡st❛❜❧✐s❤ ❛♥ ✐♠♣♦rt❛♥t ♣r♦♣❡rt② ♦❢ t❤❡ P● ( b , c ) ❞❡♥s✐t② ✕ ♥❛♠❡❧②✱ t❤❛t ❢♦r a ∈ ℜ ❛♥❞ η ∈ ℜ ✱ � ∞ ( ❡ η ) a ❡ − ωη ✷ / ✷ p ( ω | b , ✵ ) ❞ ω, ( ✶ + ❡ η ) b = ✷ − b ❡ κη ✵ ✇❤❡r❡ κ = a − b / ✷ ❛♥❞ p ( ω | b , ✵ ) ❞❡♥♦t❡s ❛ P● ( b , ✵ ) ❞❡♥s✐t②✳ ❚❤❡ ❧❡❢t✲❤❛♥❞ s✐❞❡ ❤❛s t❤❡ s❛♠❡ ❢✉♥❝t✐♦♥❛❧ ❢♦r♠ ❛s t❤❡ ♣r♦❜❛❜✐❧✐t② ♣❛r❛♠❡t❡r ✐♥ ❛ ❧♦❣✐st✐❝ r❡❣r❡ss✐♦♥ ♠♦❞❡❧ ❚❤❡ ✐♥t❡❣r❛♥❞ ♦♥ t❤❡ r✐❣❤t✲❤❛♥❞ s✐❞❡ ✐s t❤❡ ❦❡r♥❡❧ ♦❢ ❛ ♥♦r♠❛❧ ❞❡♥s✐t② ✇✐t❤ ♣r❡❝✐s✐♦♥ ω ✭✐✳❡✳✱ t❤❡ ❝♦♥❞✐t✐♦♥❛❧ ❞❡♥s✐t② ♦❢ η ✮ t✐♠❡s t❤❡ ♣r✐♦r ❢♦r ω ✷✾ ✴ ✶✻✵

❈♦♥♥❡❝t✐♦♥ t♦ t❤❡ ▲♦❣✐st✐❝ ▼♦❞❡❧ ■♥ ♣❛rt✐❝✉❧❛r✱ ✉♥❞❡r t❤❡ ❧♦❣✐st✐❝ ♠♦❞❡❧✱ t❤❡ ❧✐❦❡❧✐❤♦♦❞ ❢♦r t❤❡ ❜✐♥❛r② r❡s♣♦♥s❡ ✈❡❝t♦r ❨ = ( Y ✶ , . . . , Y n ) T ✐s � n p ( ② | β ) = p ( y i | β ) i = ✶ � � y i � � ✶ − y i � n ❡①♣ ( η i ) ✶ = ✶ + ❡①♣ ( η i ) ✶ + ❡①♣ ( η i ) i = ✶ � n ( ❡ η i ) y i = ✶ + ❡ η i , i = ✶ ✇❤❡r❡ η i = ① T i β ✳ ❚❤❡ i ✲t❤ ❡❧❡♠❡♥t ♦❢ t❤❡ ❇❡r♥♦✉❧❧✐ ❧✐❦❡❧✐❤♦♦❞ ❤❛s t❤❡ s❛♠❡ ❢♦r♠ ❛s t❤❡ ❧❡❢t✲❤❛♥❞ ❡①♣r❡ss✐♦♥ ✐♥ t❤❡ ❡❛r❧✐❡r ♣r♦♣❡rt②✱ ✇✐t❤ a i = Y i ❛♥❞ b = ✶✳ ✸✵ ✴ ✶✻✵

❈♦♥♥❡❝t✐♦♥ t♦ t❤❡ ▲♦❣✐st✐❝ ▼♦❞❡❧ ❚❤✉s✱ ✇❡ ❝❛♥ r❡✲✇r✐t❡ t❤❡ ❇❡r♥♦✉❧❧✐ ❧✐❦❡❧✐❤♦♦❞ ✐♥ t❡r♠s ♦❢ t❤❡ Pó❧②❛✲●❛♠♠❛ r❛♥❞♦♠ ✈❛r✐❛❜❧❡s ω = ( ω ✶ , . . . , ω n ) T ✿ � ∞ ❡ − ω i η ✷ p ( y i | β ) ∝ ❡ κ i η i i / ✷ p ( ω i | ✶ , ✵ ) ❞ ω i , ✵ ✇❤❡r❡ κ i = y i − ✶ ✷ ✳ ❙♦✱ ❛t ❡❛❝❤ ✐t❡r❛t✐♦♥✱ ✇❡✬❧❧ ❞r❛✇ ω i ✱ ✉♣❞❛t❡ β ❜❛s❡❞ ♦♥ t❤❡ ♥♦r♠❛❧ ♠♦❞❡❧ ❢♦r η i = ① T i β ✱ ❛♥❞ t❤❡♥ ❛✈❡r❛❣❡ ❛❝r♦ss t❤❡ ✐t❡r❛t✐♦♥s t♦ ♣❡r❢♦r♠ t❤❡ ▼♦♥t❡ ❈❛r❧♦ ✐♥t❡❣r❛t✐♦♥ ✸✶ ✴ ✶✻✵

❉✐str✐❜✉t✐♦♥ ♦❢ ▲❛t❡♥t ◆♦r♠❛❧s ❩ ❇② ❛♣♣❡❛❧✐♥❣ t♦ t❤❡ ❛❜♦✈❡ ♣r♦♣❡rt✐❡s t❤❡ Pó❧②❛✲●❛♠♠❛ ❞✐str✐❜✉t✐♦♥✱ P♦❧s♦♥ ❡t ❛❧✳ s❤♦✇ t❤❛t t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❞✐str✐❜✉t✐♦♥ ♦❢ β ✱ ❣✐✈❡♥ ❨ ❛♥❞ ω ✱ ✐s � � − ✶ ✷ ( ③ − ❳ β ) T ❲ ( ③ − ❳ β ) p ( β | ❨ = ② , ω ) ∝ π ( β ) ❡①♣ , ✇❤❡r❡ • π ( β ) ✐s t❤❡ ♣r✐♦r ❞✐str✐❜✉t✐♦♥ ❢♦r β • ❋♦r i = ✶ , . . . , n ✱ z i = y i − ✶ / ✷ ✇✐t❤ ③ = ( z ✶ , . . . , z n ) T ω i • ❲ = ❞✐❛❣ ( ω i ) ✐s ❛♥ n × n ♣r❡❝✐s✐♦♥ ♠❛tr✐① ■t ✐s ❝❧❡❛r t❤❛t t❤❡ r❛♥❞♦♠ ✈❛r✐❛❜❧❡ ❩ ✐s ♥♦r♠❛❧❧② ❞✐str✐❜✉t❡❞ ✇✐t❤ ♠❡❛♥ η = ❳ β ❛♥❞ ❞✐❛❣♦♥❛❧ ❝♦✈❛r✐❛♥❝❡ ❲ − ✶ ✸✷ ✴ ✶✻✵

❋✉❧❧ ❈♦♥❞✐t✐♦♥❛❧ ❢♦r β ❚❤✉s✱ ❛ss✉♠✐♥❣ ❛ ◆ p ( β ✵ , ❚ − ✶ ✵ ) ♣r✐♦r ❢♦r β ✱ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r β ❣✐✈❡♥ ❩ = ③ ❛♥❞ ❲ ✐s ◆ p ( ♠ , ❱ ) ✱ ✇❤❡r❡ � � − ✶ ❚ ✵ + ❳ T ❲ ❳ = ❱ � � ❚ ✵ β ✵ + ❳ T ❲ ③ ♠ = ❱ ❚❤✐s ❧❡❛❞s t♦ t❤❡ ❢♦❧❧♦✇✐♥❣ ●✐❜❜s s❛♠♣❧❡r✿ ✶ ❋♦r i = ✶ , . . . , n ✱ ✉♣❞❛t❡ ω i ❢r♦♠ ❛ P● ( ✶ , η i ) ❞❡♥s✐t②✱ ✇❤❡r❡ η i = ① T i β ✷ ❋♦r i = ✶ , . . . , n ✱ ❞❡✜♥❡ z i = y i − ✶ / ✷ ω i ✸ ❈♦♥❞✐t✐♦♥❛❧ ♦♥ ③ ❛♥❞ ❲ ✱ ✉♣❞❛t❡ β ❢r♦♠ ◆ p ( ♠ , ❱ ) ✱ ✇❤❡r❡ ♠ ❛♥❞ ❱ ❛r❡ ❣✐✈❡♥ ❛❜♦✈❡✳ ✸✸ ✴ ✶✻✵

❙❛♠♣❧✐♥❣ ❢r♦♠ t❤❡ P● ❉❡♥s✐t② ❆❝❝❡♣t❛♥❝❡ r❡❥❡❝t✐♦♥ s❛♠♣❧✐♥❣ ✐s ✉s❡❞ t♦ ❞r❛✇ ❢r♦♠ t❤❡ P● ❞❡♥s✐t② ❚❤✐s ❝❛♥ ❜❡ ✐♠♣❧❡♠❡♥t❡❞ ✉s✐♥❣ t❤❡ r♣❣ ❢✉♥❝t✐♦♥ ❢r♦♠ t❤❡ ❘ ♣❛❝❦❛❣❡ ❇❛②❡s▲♦❣✐t ■ ❤❛✈❡ ✉♣❧♦❛❞❡❞ ❛ ③✐♣ ✜❧❡ ♦❢ t❤❡ ♣❛❝❦❛❣❡ ♦♥t♦ t❤❡ ❝♦✉rs❡ ✇❡❜s✐t❡ ❆❧t❡r♥❛t✐✈❡❧②✱ ②♦✉ ❝❛♥ ❞♦✇♥❧♦❛❞ ❢r♦♠ ❤tt♣s✿✴✴♠r❛♥✳♠✐❝r♦s♦❢t✳❝♦♠✴s♥❛♣s❤♦t✴✷✵✶✹✲✶✵✲✷✵✴ ✇❡❜✴♣❛❝❦❛❣❡s✴❇❛②❡s▲♦❣✐t✴✐♥❞❡①✳❤t♠❧ ❖r ❣♦ t♦ ❤tt♣s✿✴✴❣✐t❤✉❜✳❝♦♠✴❥✇✐♥❞❧❡✴❇❛②❡s▲♦❣✐t ❢♦r ❛ ♠♦r❡ r❡❝❡♥t ✈❡rs✐♦♥ ❙❡❡ P● ▲♦❣✐t✳r ❢♦r ❞❡t❛✐❧s ✸✹ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r Pó❧②❛✲●❛♠♠❛ ▲♦❣✐st✐❝ ▼♦❞❡❧ Logistic Regression Using PG Latent Variables library (BayesLogit) # For rpg function library (mvtnorm) # Priors beta0<- rep (0,p) # Prior mean of beta T0<- diag (.01,p) # Prior precision of beta # Inits beta<- rep (0,p) ################# # Store Samples # ################# nsim<-1000 # Number of MCMC Iterations thin<-1 # Thinning interval burn<-nsim/2 # Burnin lastit<-(nsim-burn)/thin # Last stored value Beta<- matrix (0,lastit,p) ######### # Gibbs # ######### tmp<- proc.time () # Store current time for (i in 1:nsim){ eta<-X%*%beta w<- rpg (n,1,eta) z<-(y-1/2)/w # Or define z=y-1/2 and omit w in posterior mean m below v<- solve ( crossprod (X* sqrt (w))+T0) # Or solve(X%*%W%*%X), where W=diag(w) -- but this is slower m<-v%*%(T0%*%beta0+ t (w*X)%*%z) # Can omit w here if you define z=y-1/2 beta<- c ( rmvnorm (1,m,v)) ################# # Store Results # ################# if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta } if (i%%100==0) print (i) } proc.time ()-tmp # MCMC run time -- 1.2 seconds to run 1000 iterations with n=1000 ✸✺ ✴ ✶✻✵

❊①❛♠♣❧❡ ❚❤❡ ♣r♦❣r❛♠ P● ▲♦❣✐t✳r ✜ts t❤❡ s❛♠❡ ❧♦❣✐st✐❝ ♠♦❞❡❧ ❛s ❜❡❢♦r❡✿ Y i ∼ ❇❡r♥ ( π i ) ❧♦❣✐t ( π i ) = β ✵ + β ✶ x i , i = ✶ , . . . , ✶✵✵✵ . ❚❤❡ r❡s✉❧ts✱ ❜❛s❡❞ ♦♥ ✶✵✵✵ ✐t❡r❛t✐♦♥s ✇✐t❤ ❛ ❜✉r♥✲✐♥ ♦❢ ✺✵✵✱ ❛r❡ ✐❞❡♥t✐❝❛❧ t♦ t❤❡ ♣r❡✈✐♦✉s ♦♥❡s✿ ❚❛❜❧❡ ✹✿ ❘❡s✉❧ts ❢♦r ▲♦❣✐st✐❝ ▼♦❞❡❧✳ ❚r✉❡ P♦st❡r✐♦r ▼❡❛♥ ✭❙❉✮ † P❛r❛♠❡t❡r ❱❛❧✉❡ ▼▲❊ ✭❙❊✮ β ✵ − ✶ − ✶ . ✵✽ ( ✵ . ✵✽ ) − ✶ . ✵✽ ( ✵ . ✵✽ ) β ✶ ✶ ✵ . ✾✷ ( ✵ . ✵✾ ) ✵ . ✾✸ ( ✵ . ✵✾ ) † ❇❛s❡❞ ♦♥ P● s❛♠♣❧❡r✳ ✸✻ ✴ ✶✻✵

❇❛②❡s✐❛♥ ❖r❞✐♥❛❧ ❘❡❣r❡ss✐♦♥ ❆❧❜❡rt ❛♥❞ ❈❤✐❜ ✭✶✾✾✸✮ ❡①t❡♥❞ t❤❡✐r ❞❛t❛ ❛✉❣♠❡♥t❛t✐♦♥ s❛♠♣❧❡r t♦ ❛❝❝♦♠♠♦❞❛t❡ ♦r❞❡r❡❞ ❝❛t❡❣♦r✐❝❛❧ ♦✉t❝♦♠❡s ❋♦r ❡①❛♠♣❧❡✱ t❤❡ ❝✉♠✉❧❛t✐✈❡ ❧♦❣✐t ♠♦❞❡❧ ✇✐t❤ K ❝❛t❡❣♦r✐❡s ✐s ❧♦❣✐t ( φ ik ) = ❧♦❣✐t [ Pr ( Y i ≤ k )] = γ k + ① T i β , k = ✶ , . . . , K − ✶ , • φ ik ✐s t❤❡ k ✲t❤ ❝✉♠✉❧❛t✐✈❡ ♣r♦❜❛❜✐❧✐t② • γ k ✐s t❤❡ ✐♥t❡r❝❡♣t ❛ss♦❝✐❛t❡❞ ✇✐t❤ k ✲t❤ ❝✉♠✉❧❛t✐✈❡ ❧♦❣✐t • ① i ✐s ❛ ✈❡❝t♦r ♦❢ ❝♦✈❛r✐❛t❡s ❡①❝❧✉❞✐♥❣ ❛♥ ✐♥t❡r❝❡♣t • β ✐s ❛ ✈❡❝t♦r ♦❢ r❡❣r❡ss✐♦♥ ❝♦❡✣❝✐❡♥ts ❝♦♠♠♦♥ t♦ ❛❧❧ ❝✉♠✉❧❛t✐✈❡ ❧♦❣✐ts ✭♣r♦♣♦rt✐♦♥❛❧ ♦❞❞s✮ β j ❞❡♥♦t❡s t❤❡ ✐♥❝r❡❛s❡ ✐♥ t❤❡ ❧♦❣ ♦❞❞s ♦❢ ❧♦✇❡r ❝❛t❡❣♦r② ✈❛❧✉❡s ❢♦r ❛ ✉♥✐t ✐♥❝r❡❛s❡ ✐♥ x j ✸✼ ✴ ✶✻✵

▲❛t❡♥t ❱❛r✐❛❜❧❡ ■♥t❡r♣r❡t❛t✐♦♥ ♦❢ ❖r❞✐♥❛❧ ▼♦❞❡❧ 0.4 0.3 Pr( γ 1 < Z i < γ 2 ) = Pr(Y i = 2) f(z i ) 0.2 0.1 Pr(Z i < γ 1 ) = Pr(Y i = 1) Pr(Z i > γ 2 ) = Pr(Y i = 3) 0.0 γ 1 γ 2 T β x i z i ✸✽ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r Pr♦♣♦rt✐♦♥❛❧ ❖❞❞s ▼♦❞❡❧ ❚❤❡ ●✐❜❜s s❛♠♣❧❡r ✐s ❛ s❧✐❣❤t ♠♦❞✐✜❝❛t✐♦♥ ♦❢ t❤❡ ♣r❡✈✐♦✉s s❛♠♣❧❡r ❢♦r ♣r♦❜✐t✴❧♦❣✐t r❡❣r❡ss✐♦♥ ❍❡r❡✱ ✐♥ ❛❞❞✐t✐♦♥ t♦ ✉♣❞❛t✐♥❣ Z i ❛♥❞ β ✱ ✇❡ ♠✉st ❛❧s♦ ✉♣❞❛t❡ γ = ( γ ✶ , . . . , γ K − ✶ ) T ◆♦t❡ t❤❛t ✇❡ ❤❛✈❡ t❤❡ ❢♦❧❧♦✇✐♥❣ ♦r❞❡r ❝♦♥str❛✐♥t✿ ❛❧❧ Z ( y = k ) < γ k < ❛❧❧ Z ( y =( k + ✶ )) ❢♦r k = ✶ , . . . , K − ✶ ❚❤✐s ✐♠♣❧✐❡s✱ ❢♦r ✐♥st❛♥❝❡✱ t❤❛t γ ✶ ❤❛s t♦ ❧✐❡ ❜❡t✇❡❡♥ t❤❡ ❧❛r❣❡st Z s✳t✳ Y = ✶ ❛♥❞ t❤❡ s♠❛❧❧❡st Z s✳t✳ t❤❛t Y = ✷ ❲❡ ❤❛✈❡ t♦ ♦❜❡② t❤✐s ❝♦♥str❛✐♥t ✇❤❡♥ ✉♣❞❛t✐♥❣ γ ✸✾ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r Pr♦♣♦rt✐♦♥❛❧ ❖❞❞s ▼♦❞❡❧ ❚❤❡ s❛♠♣❧❡r ♣r♦❝❡❡❞s ❜② ✐♥✐t✐❛❧✐③✐♥❣ γ ❛♥❞ β ❛♥❞ t❤❡♥ ✐t❡r❛t✐♥❣ t❤r♦✉❣❤ t❤❡ ❢♦❧❧♦✇✐♥❣ st❡♣s✿ i β , λ − ✶ ✶ ❋♦r ✐❂✶✱✳ ✳ ✳ ✱♥✱ ✐❢ y i = k ✱ ❞r❛✇ z i ❢r♦♠ ◆ ( ① T i ) tr✉♥❝❛t❡❞ ❜❡t✇❡❡♥ γ k − ✶ ❛♥❞ γ k ✱ ✇✐t❤ γ ✵ = −∞ ❛♥❞ γ K = ∞ ◆♦t❡✿ ❋♦r t❤❡ ❝✉♠✉❧❛t✐✈❡ ♣r♦❜✐t ♠♦❞❡❧✱ λ i = ✶ ∀ i ✷ ❯♣❞❛t❡ β ❢r♦♠ ◆ p ( ♠ , ❱ ) ✱ � � ❚ ✵ + ❳ T ❲ ❳ ❱ = � � ❚ ✵ β ✵ + ❳ T ❲ ③ ♠ = ❱ ✇❤❡r❡ ❲ = ■ n ❢♦r ♣r♦❜✐t ♠♦❞❡❧ ❛♥❞ ❲ = ❞✐❛❣ ( λ i ) ❢♦r t ✲❧✐♥❦ ♠♦❞❡❧ ✸ ❢♦r k = ✶ , . . . , K − ✶✱ ✉♣❞❛t❡ γ k ❢r♦♠ ❯♥✐❢ ( ♠❛① z i : y i = k , ♠✐♥ z i : y i = k + ✶ ) ✹ ❋♦r t❤❡ ❝✉♠✉❧❛t✐✈❡ ❧♦❣✐t ♠♦❞❡❧✱ ✉♣❞❛t❡ λ i ❢r♦♠ ✐ts ●❛♠♠❛ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❛♥❛❧♦❣♦✉s t♦ t❤❡ ✭✷✲❝❛t❡❣♦r②✮ ❧♦❣✐st✐❝ ♠♦❞❡❧ ❞✐s❝✉ss❡❞ ❡❛r❧✐❡r ❙❡❡ ❈✉♠❴▲♦❣✐t✳r ❛♥❞ ❈✉♠❴Pr♦❜✐t✳r ❢♦r ❞❡t❛✐❧s ✹✵ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r ❈✉♠✉❧❛t✐✈❡ ▲♦❣✐t ▼♦❞❡❧ Gibbs Sampler for 3-Category Cumulative Logit Model First assign a N ( β 0 , T − 1 0 ) prior to β # Initial Values gam1<-0 gam2<-3 beta<-0 # Only 1 covariate in this example lambda<- rep (1,n) # Weights nu<-8 # 8 df ~ logistic -- assume fixed z<- rep (0,n) # latent z vector ################### # GIBBS SAMPLER # ################### tmp<- proc.time () for (i in 1:nsim) { # Draw latent z using inverse CDF method muz<-x*beta z[y==1]<- qnorm ( runif (ns[1],0, pnorm (gam1,muz[y==1], sqrt (1/lambda[y==1]))),muz[y==1], sqrt (1/lambda[y==1 z[y==2]<- qnorm ( runif (ns[2], pnorm (gam1,muz[y==2], sqrt (1/lambda[y==2])), pnorm (gam2,muz[y==2], sqrt (1/lambda[y== z[y==3]<- qnorm ( runif (ns[3], pnorm (gam2,muz[y==3], sqrt (1/lambda[y==3])),1),muz[y==3], sqrt (1/lambda[y==3 # Update gammas gam1<- runif (1, max (z[y==1]), min (z[y==2])) gam2<- runif (1, max (z[y==2]), min (z[y==3])) # Update beta vbeta<- solve ( crossprod (x* sqrt (lambda))) mbeta<-vbeta%*%( crossprod (x*lambda,z)) beta<- rnorm (1,mbeta, sqrt (vbeta)) # Update lambda lambda<- rgamma (n,(nu+1)/2,(nu+(z-x*beta)^2)/2) ################# # Store Results # ################# if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j]<-beta Gam1[j]<-gam1 Gam2[j]<-gam2 } if(i%%100==0) print (i) } proc.time ()-tmp # MCMC run time -- took 4.5 seconds to run 10K iterations with n=500 ✹✶ ✴ ✶✻✵

❊①❛♠♣❧❡ ❚❤❡ ♣r♦❣r❛♠ ❈✉♠❴▲♦❣✐t✳r s✐♠✉❧❛t❡s ✺✵✵ ♦❜s❡r✈❛t✐♦♥s ❛❝❝♦r❞✐♥❣ t♦ t❤❡ t❤❡ ❢♦❧❧♦✇✐♥❣ t❤r❡❡✲❝❛t❡❣♦r② ❝✉♠✉❧❛t✐✈❡ ❧♦❣✐t ♠♦❞❡❧✿ ❧♦❣✐t ( φ ik ) = ❧♦❣✐t [ Pr ( Y i ≤ k )] = γ k + β x i , k = ✶ , ✷ . ❚❤❡ r❡s✉❧ts✱ ❜❛s❡❞ ♦♥ ✶✵✱✵✵✵ ✐t❡r❛t✐♦♥s ✇✐t❤ ❛ ❜✉r♥✲✐♥ ♦❢ ✺✵✵✵ ❛♥❞ t❤✐♥♥✐♥❣ ♦❢ ✶✵✱ ❛r❡✿ ❚❛❜❧❡ ✺✿ ❘❡s✉❧ts ❢♦r ❈✉♠✉❧❛t✐✈❡ ▲♦❣✐t ▼♦❞❡❧✳ ❚r✉❡ P♦st❡r✐♦r ▼❡❛♥ ✭❙❉✮ † P❛r❛♠❡t❡r ❱❛❧✉❡ ▼▲❊ ✭❙❊✮ γ ✶ − ✶ − ✶ . ✵✹ ( ✵ . ✶✼ ) − ✶ . ✵✵ ( ✵ . ✶✾ ) γ ✷ ✶ ✵ . ✽✼ ( ✵ . ✶✻ ) ✵ . ✾✺ ( ✵ . ✶✼ ) β ✵ . ✼✺ ✵ . ✼✸ ( ✵ . ✵✽ ) ✵ . ✼✼ ( ✵ . ✵✽ ) † ❇❛s❡❞ ♦♥ ❆❧❜❡rt ❛♥❞ ❈❤✐❜ s❛♠♣❧❡r✳ ✹✷ ✴ ✶✻✵

❚r❛❝❡ P❧♦ts ❋✐❣✉r❡ ✶✿ ❚r❛❝❡ ♣❧♦ts s❤♦✇ ❤✐❣❤ ❛✉t♦❝♦rr❡❧❛t✐♦♥✱ ❡s♣✳ ❛♠♦♥❣ γ s −0.6 −1.0 γ 1 −1.4 0 100 200 300 400 500 Iteration 1.2 1.0 γ 2 0.8 0.6 0 100 200 300 400 500 Iteration 0.8 β 0.6 0 100 200 300 400 500 Iteration ✹✸ ✴ ✶✻✵

▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ❘❡❣r❡ss✐♦♥ ▲❡t Y i ❞❡♥♦t❡ ❛♥ ✉♥♦r❞❡r❡❞ ♦r ♥♦♠✐♥❛❧ ❝❛t❡❣♦r✐❝❛❧ ❘❱ t❛❦✐♥❣ ✈❛❧✉❡s k = ✶ , . . . , K ✇✐t❤ k ✲t❤ ♣r♦❜❛❜✐❧✐t② π ik = Pr ( Y i = k | ① i ) ✱ ✇❤❡r❡ � K k = ✶ π ik = ✶ ❚❤❡ ♠✉❧t✐♥♦♠✐❛❧ ❧♦❣✐t ♠♦❞❡❧ ✐s ❣✐✈❡♥ ❜② ❡ ① T i β k Pr ( Y i = k | ① i ) = π ik = i β j , k = ✶ , . . . , K , � K j = ✶ ❡ ① T ✇✐t❤ β K = ✵ ❢♦r t❤❡ r❡❢❡r❡♥❝❡ ❝❛t❡❣♦r② K ✳ ❆❧t❡r♥❛t✐✈❡❧②✱ ✇❡ ❝❛♥ ✇r✐t❡ ❡ ① T i β k Pr ( Y i = k | ① i ) = π ik = i β j , k = ✶ , . . . , K − ✶ ✶ + � K − ✶ j = ✶ ❡ ① T ✶ Pr ( Y i = K | ① i ) = π iK = ✶ + � K − ✶ r❡❢❡r❡♥❝❡ ❝❛t❡❣♦r② , j = ✶ ❡ ① T i β j ✇❤❡r❡ β k ❞❡♥♦t❡s t❤❡ r❡❣r❡ss✐♦♥ ❝♦❡✣❝✐❡♥ts ❛ss♦❝✐❛t❡❞ ✇✐t❤ k ✲t❤ ❝❛t❡❣♦r② ■♥ ♦t❤❡r ✇♦r❞s✱ ❝♦♥❞✐t✐♦♥❛❧ ♦♥ ① i ✱ Y i ❤❛s ❛ ▼✉❧t✐ ( ✶ , π i ) ✱ ✇❤❡r❡ π i = ( π i ✶ , . . . , π iK ) ❛♥❞ � K k = ✶ π ik = ✶✳ ✹✹ ✴ ✶✻✵

▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ❘❡❣r❡ss✐♦♥ ❲❡ ❝❛♥ ❛❧s♦ ✇r✐t❡ t❤❡ ♠✉❧t✐♥♦♠✐❛❧ ❧♦❣✐t ♠♦❞❡❧ ❛s Pr ( Y i = K | ① i ) = π ik Pr ( Y i = k | ① i ) = ❡①♣ ( ① T i β k ) , k = ✶ , . . . , K − ✶ , π iK ✇❤✐❝❤ ✐s t❤❡ ♦❞❞s ♦❢ ❜❡✐♥❣ ✐♥ ❝❛t❡❣♦r② k ✈s✳ t❤❡ r❡❢❡r❡♥❝❡ ❝❛t❡❣♦r② K ✳ ❖r✱ ❡q✉✐✈❛❧❡♥t❧②✱ � π ik � = ① T ❧♥ i β k . π iK ■♥ t❤❡ ♣r♦♣♦rt✐♦♥❛❧ ♦❞❞s ♠♦❞❡❧✱ ♦♥❧② t❤❡ ✐♥t❡r❝❡♣ts α k ✈❛r✐❡❞ ❛❝r♦ss ❝❛t❡❣♦r✐❡s✱ ✇❤✐❧❡ β r❡♠❛✐♥❡❞ ❝♦♥st❛♥t ❍❡r❡✱ ❛❧❧ r❡❣r❡ss✐♦♥ ♣❛r❛♠❡t❡rs ✈❛r② ✇✐t❤ r❡s♣❡❝t t♦ k ✭❤❡♥❝❡ β k ✮ ❚❤❡r❡ ❛r❡ ❛❧s♦ ♣❛rt✐❛❧ ♣r♦♣♦rt✐♦♥❛❧ ♦❞❞s ♠♦❞❡❧s t❤❛t ❛❧❧♦✇ ♦♥❧② s♦♠❡ β ✬s t♦ ✈❛r② ❛❝r♦ss ❝❛t❡❣♦r✐❡s ✹✺ ✴ ✶✻✵

▲❛t❡♥t ❱❛r✐❛❜❧❡ ■♥t❡r♣r❡t❛t✐♦♥s ❇♦t❤ t❤❡ ♣r♦♣♦rt✐♦♥❛❧ ♦❞❞s ❛♥❞ ♠✉❧t✐♥♦♠✐❛❧ ❧♦❣✐t ♠♦❞❡❧s ❤❛✈❡ ❧❛t❡♥t ✈❛r✐❛❜❧❡ ✐♥t❡r♣r❡t❛t✐♦♥s ❋♦r t❤❡ ▼◆▲ ♠♦❞❡❧✱ ✇❡ ❤❛✈❡ K ✉♥❞❡r❧②✐♥❣ ✏✉t✐❧✐t✐❡s✑✿ = X i β ✶ + e i ✶ U i ✶ ✳ ✳ ✳ ✳ ✳ ✳ = X i β ( K − ✶ ) + e i ( K − ✶ ) U i ( K − ✶ ) U iK = e iK , ✇❤❡r❡ t❤❡ e ik ✬s ❢♦❧❧♦✇ ✐✳✐✳❞✳ st❛♥❞❛r❞ ●✉♠❜❡❧ ♦r ❚②♣❡✲■ ❡①tr❡♠❡ ✈❛❧✉❡ ❞✐str✐❜✉t✐♦♥s ✇✐t❤ ❈❉❋s ♦❢ t❤❡ ❢♦r♠ ❋ X ( x ) = ❡ − ❡ − x ❙❡t Y i = k ✐✳❢✳❢✳ U ik = ♠❛① ( U i ✶ , . . . , U iK ) ❍❲✿ ❙❤♦✇ t❤❛t t❤❡ ❞✐✛❡r❡♥❝❡ ♦❢ t✇♦ ✐♥❞❡♣❡♥❞❡♥t st❛♥❞❛r❞ ●✉♠❜❡❧ ❘✳❱✳s ❢♦❧❧♦✇s ❛ ❧♦❣✐st✐❝ ❞✐str✐❜✉t✐♦♥ ✹✻ ✴ ✶✻✵

❉✐s❝r❡t❡ ❈❤♦✐❝❡ ▼♦❞❡❧s ❚❤✐s ❢♦r♠s t❤❡ ❜❛s✐s ♦❢ s♦✲❝❛❧❧❡❞ r❛♥❞♦♠ ✉t✐❧✐t② ♦r ❞✐s❝r❡t❡ ❝❤♦✐❝❡ ♠♦❞❡❧s ❙❡❡ ❚r❛✐♥ ❉✐s❝r❡t❡ ❈❤♦✐❝❡ ▼❡t❤♦❞s ✇✐t❤ ❙✐♠✉❧❛t✐♦♥ ❢♦r ♠♦r❡ ❞❡t❛✐❧s ♦♥ ❞✐s❝r❡t❡ ❝❤♦✐❝❡ ♠♦❞❡❧s ❚❤❡ ❛✉t❤♦r ❞✐s❝✉ss❡s ♣r♦s ❛♥❞ ❝♦♥s ♦❢ t❤❡ ♠✉❧t✐♥♦♠✐❛❧ ❧♦❣✐t ♠♦❞❡❧✱ ✐♥❝❧✉❞✐♥❣ t❤❡ ✐♥❞❡♣❡♥❞❡♥❝❡ ♦❢ ✐rr❡❧❡✈❛♥t ❛❧t❡r♥❛t✐✈❡s ✭■■❆✮ ♣r♦♣❡rt② ❛s ❡①❡♠♣❧✐✜❡❞ ✐♥ t❤❡ ❝❧❛ss✐❝ ❘❡❞ ❇✉s✴❇❧✉❡ ❇✉s ♣r♦❜❧❡♠ ❚♦ ❝✐r❝✉♠✈❡♥t t❤✐s ♣r♦❜❧❡♠✱ ♦♥❡ ❝❛♥ ✉s❡ t❤❡ ♥❡st❡❞ ❧♦❣✐t ♦r ♠✉❧t✐♥♦♠✐❛❧ ♣r♦❜✐t ♠♦❞❡❧ ❋♦r ❢✉rt❤❡r ❞❡t❛✐❧s✱ s❡❡ ❆❣r❡st✐ ❈❛t❡❣♦r✐❝❛❧ ❉❛t❛ ❆♥❛②s✐s ❙❡❝t✐♦♥s ✽✳✺ ❛♥❞ ✽✳✻ ❛♥❞ r❡❢❡r❡♥❝❡s t❤❡r❡✐♥ ✹✼ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ▼♦❞❡❧ ❚❤❡ Pó❧②❛✲●❛♠♠❛ ❞❛t❛✲❛✉❣♠❡♥t❛t✐♦♥ ❛♣♣r♦❛❝❤ ❝❛♥ ❜❡ ❡①t❡♥❞❡❞ t♦ ❤❛♥❞❧❡ ♠✉❧t✐♥♦♠✐❛❧ ❧♦❣✐t r❡❣r❡ss✐♦♥ ❘❡❝❛❧❧ t❤❛t t❤❡ ♠✉❧t✐♥♦♠✐❛❧ ❧♦❣✐t ♠♦❞❡❧ ❝❛♥ ❜❡ ✇r✐tt❡♥ ❛s ❡ ① T i β k Pr ( Y i = k | ① i ) = π ik = i β j , k = ✶ , . . . , K , � K j = ✶ ❡ ① T ✇✐t❤ β K = ✵ ❢♦r t❤❡ r❡❢❡r❡♥❝❡ ❝❛t❡❣♦r② K ✳ ▲❡t✬s ❝♦♥s✐❞❡r t❤❡ ✉♣❞❛t❡ ❢♦r β k ( k = ✶ , . . . , K − ✶ ) ❝♦♥❞✐t✐♦♥❛❧ ♦♥ ❜♦t❤ ② ❛♥❞ β j ∀ j � = k ✳ ❚❤❡ ✐❞❡❛ ✐s t♦ ❝②❝❧❡ t❤r♦✉❣❤ t❤❡ ✉♣❞❛t❡s ♦❢ β k ♦♥❡ ❛t ❛ t✐♠❡ ❝♦♥❞✐t✐♦♥❛❧ ♦♥ t❤❡ ♦t❤❡r β ✬s ✹✽ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ▼♦❞❡❧ ∗ ❋♦❧❧♦✇✐♥❣ ❍♦❧♠❡s ❛♥❞ ❍❡❧❞ ✭✷✵✵✻✮✱ ✇❡ ❝❛♥ ✇r✐t❡ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r β k ✱ ❣✐✈❡♥ ② ❛♥❞ β j � = k ✱ ❛s ❛ ❢✉♥❝t✐♦♥ ♦❢ ❛ ❇❡r♥♦✉❧❧✐ ❧✐❦❡❧✐❤♦♦❞ � n π U ik ik ( ✶ − π ik ) ✶ − U ik f ( β k | ② , β j � = k ) ∝ f ( β k ) i = ✶ ✇❤❡r❡ • f ( β k ) ❞❡♥♦t❡s t❤❡ ♣r✐♦r ❢♦r β k • U ik = ✶ ( Y i = k ) ✐s ❛♥ ✐♥❞✐❝❛t♦r t❤❛t Y i = k ❡ ① T i β k • π ik = Pr ( Y i = k ) = Pr ( U ik = ✶ ) = j = ✶ ❡ ① T i β j � K ✹✾ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ▼♦❞❡❧ ∗ ❋✉rt❤❡r✱ ✇❡ ❝❛♥ r❡✇r✐t❡ π ik ❛s ❡ ① T i β k − c ik π ik = Pr ( U ik = ✶ ) = ✶ + ❡ ① T i β k − c ik ❡ η ik = ✶ + ❡ η ik , ✇❤❡r❡ c ik = ❧♦❣ � j � = k ❡ ① T i β j ❛♥❞ η ik = ① T i β k − c ik ✳ ◆♦t❡ t❤❛t � j � = k ❡ ① T i β j ✐♥❝❧✉❞❡s t❤❡ r❡❢❡r❡♥❝❡ ❝❛t❡❣♦r② K ✳ ❇❡❝❛✉s❡ β K = ✵ ✱ ✐t ❢♦❧❧♦✇s t❤❛t ❡ ① T i β K = ✶ ❛♥❞ ❤❡♥❝❡   � � ❡ ① T i β j = ❧♦❣  ✶ + ❡ ① T  i β j c ik = ❧♦❣ j � = k j / ∈{ k , K } ✺✵ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ▼♦❞❡❧ ❚❤✉s✱ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r β k ❣✐✈❡♥ ② ❛♥❞ β j � = k ✱ ✐s � � U ik � � ✶ − U ik � n ❡ η ik ✶ f ( β k | ② , β j � = k ) ∝ f ( β k ) ✶ + ❡ η ik ✶ + ❡ η ik i = ✶ � n ( ❡ η ik ) U ik = f ( β k ) ✶ + ❡ η ik , i = ✶ ✇❤✐❝❤ ✐s ❡ss❡♥t✐❛❧❧② ❛ ❧♦❣✐st✐❝ r❡❣r❡ss✐♦♥ ❧✐❦❡❧✐❤♦♦❞✳ ❚❤✉s✱ ✇❡ ❝❛♥ ❛♣♣❧② t❤❡ Pó❧②❛✲●❛♠♠❛ ❞❛t❛✲❛✉❣♠❡♥t❛t✐♦♥ s❝❤❡♠❡ t♦ ✉♣❞❛t❡ t❤❡ β k ✬s ( k = ✶ , . . . , K − ✶ ) ♦♥❡ ❛t ❛ t✐♠❡ ❜❛s❡❞ ♦♥ t❤❡ ❜✐♥❛r② ✐♥❞✐❝❛t♦rs U ik = ✶ ( Y i = k ) ✺✶ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ▼♦❞❡❧ ❆ss✉♠✐♥❣ ❛ ◆ ( β ✵ , ❚ − ✶ ✵ ) ❢♦r β ✶ , . . . , β K − ✶ ✱ t❤❡ ●✐❜❜s s❛♠♣❧❡r ♣r♦❝❡❡❞s ❛s ❢♦❧❧♦✇s✿ ✶ ❖✉ts✐❞❡ ♦❢ t❤❡ ●✐❜❜s ❧♦♦♣✱ ❞❡✜♥❡ u ik = ✶ ( y i = k ) ❢♦r i = ✶ , . . . , n ❛♥❞ k = ✶ , . . . , K − ✶ ✷ ❋♦r i = ✶ , . . . , n ❛♥❞ k = ✶ , . . . , K − ✶✱ ✉♣❞❛t❡ ω ik ❢r♦♠ ❛ P● ( ✶ , η ik ) ❞❡♥s✐t②✱ ✇❤❡r❡ η ik = ① T i β k − c ik ✱ ❛♥❞ c ik ✇❛s ❞❡✜♥❡❞ ❡❛r❧✐❡r ✸ ❋♦r i = ✶ , . . . , n ❛♥❞ k = ✶ , . . . , K − ✶✱ ❞❡✜♥❡ z ik = u ik − ✶ / ✷ + c ik ω ik ❛♥❞ ❧❡t ③ k = ( z ✶ k , . . . , z nk ) T ✹ ❋♦r k = ✶ , . . . , K − ✶✱ ✉♣❞❛t❡ β k ❢r♦♠ ◆ p ( ♠ k , ❱ k ) ✱ ✇❤❡r❡ � � − ✶ ❚ ✵ + ❳ T ❲ k ❳ ❱ k = � � ❚ ✵ β ✵ + ❳ T ❲ k ③ k ♠ k = ❱ k , ❛♥❞ ❲ k = ❞✐❛❣ ( ω ik ) ✳ ❙❡❡ P♦❧s♦♥ ❡t ❛❧✳ ✭✷✵✶✷✮ ♣❛❣❡ ✹✶ ❢♦r ❞❡t❛✐❧s✱ ❜✉t ♥♦t❡ t②♣♦ ✐♥ ❡①♣r❡ss✐♦♥ ❢♦r ♠ j ❛t ❜♦tt♦♠ ✺✷ ✴ ✶✻✵

❊①❛♠♣❧❡ ❚❤❡ ♣r♦❣r❛♠ ♠✉❧t✐♥♦♠✐❛❧✳r s✐♠✉❧❛t❡s ✶✵✵✵ ♦❜s❡r✈❛t✐♦♥s ❢r♦♠ t❤❡ ❢♦❧❧♦✇✐♥❣ t❤r❡❡✲❝❛t❡❣♦r② ♠✉❧t✐♥♦♠✐❛❧ ♠♦❞❡❧ ❡ ① T i β k Pr ( Y i = k | ① i ) = π ik = i β j , k = ✶ , . . . , ✸ , � K j = ✶ ❡ ① T ✇✐t❤ β ✶ = ✵ ❢♦r r❡❢❡r❡♥❝❡ ❝❛t❡❣♦r② ✶✳ ❚❤❡ ♣r♦❣r❛♠ ❣❡♥❡r❛t❡s ❞❛t❛ ✉s✐♥❣ t❤❡ ●✉♠❜❡❧ ❧❛t❡♥t ✉t✐❧✐t✐❡s ❙❡❡ ♠✉❧t✐♥♦♠✐❛❧✳r ❢♦r ❞❡t❛✐❧s ✺✸ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r ●❡♥❡r❛t✐♥❣ ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ❉❛t❛ Data Generation for Multinomial Logit Model # Multinomial.r # Generate and fit a 3-Category Multinomial Logit # Generate data using latent Gumbel random utilities # See Polson et al. 2012 Appendix S6.3 BUT NOTE TYPO AT BOTTOM OF PAGE 41! # 3-Category outcome: Independent, Republican and Democrat, with Ind as ref group ####################################### library (QRM) # For rGumbel function library (nnet) # To fit multinom function library (BayesLogit) # For rpg function ############################## # Generate Data under # # Random Utility Model (RUM) # ############################## set.seed (060817) n<-1000 K<-3 # Number of response categories female<- rbinom (n,1,.5) X<- cbind (1,female) p<- ncol (X) beta2<- c (1,-.5) # Males much more likely to be Rep than Ind and females "somewhat" more likely beta3<- c (.5,.5) # Males somewhat more likely to be Dem than Indep and females much more likely eta2<-X%*%beta2 eta3<-X%*%beta3 u1<- rGumbel (n, mu = 0, sigma = 1) u2<- rGumbel (n, mu = eta2, sigma = 1) u3<- rGumbel (n, mu = eta3, sigma = 1) U<- cbind (u1,u2,u3) y<- c ( apply (U,1,which.max)) # Y=1 if Ind (Reference), 2 if Rep, 3 if Dem fit<- multinom (y~female) ✺✹ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r ✸✲❈❛t❡❣♦r② ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ▼♦❞❡❧ Gibbs Sampler for Multinomial Logit Model # Define category-specific binary responses (Note: cat 1 is reference) u2<-1*(y==2) u3<-1*(y==3) # Priors beta0<- rep (0,p) # Prior mean of beta T0<- diag (.01,p) # Prior precision of beta # Inits beta2<-beta3<- rep (0,p) ################# # Gibbs sampler # ################# tmp<- proc.time () # Store current time for (i in 1:nsim){ # Update category 2 c2<- log (1+ exp (X%*%beta3)) # Note that for Cat 1, beta1=0 so that exp(X%*%beta1)= 1 eta2<-X%*%beta2-c2 w2<- rpg (n,1,eta2) z2<-(u2-1/2)/w2+c2 # Note plus sign before c2 -- Polson has typo # Could also define z2<-u2-1/2+w2*c2 and omit w2 in post. mean m v<- solve (T0+ crossprod (X* sqrt (w2))) m<-v%*%(T0%*%beta0+ t (w2*X)%*%z2) beta2<- c ( rmvnorm (1,m,v)) # Update category 3 c3<- log (1+ exp (X%*%beta2)) eta3<-X%*%beta3-c3 w3<- rpg (n,1,eta3) z3<-(u3-1/2)/w3+c3 v<- solve (T0+ crossprod (X* sqrt (w3))) m<-v%*%(T0%*%beta0+ t (w3*X)%*%z3) beta3<- c ( rmvnorm (1,m,v)) # Store if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<- c (beta2,beta3) } if (i%%100==0) print (i) } proc.time ()-tmp # MCMC run time = 2.5 seconds to run 1000 iterations ✺✺ ✴ ✶✻✵

❘❡s✉❧ts ❚❤❡ r❡s✉❧ts✱ ❜❛s❡❞ ♦♥ ✶✵✵✵ ✐t❡r❛t✐♦♥s ✇✐t❤ ❛ ❜✉r♥✲✐♥ ♦❢ ✺✵✵✱ ❛r❡✿ ❚❛❜❧❡ ✻✿ ❘❡s✉❧ts ❢♦r ▼✉❧t✐♥♦♠✐❛❧ ▲♦❣✐t ▼♦❞❡❧✳ ❚r✉❡ P♦st❡r✐♦r ▼❡❛♥ ✭❙❉✮ † P❛r❛♠❡t❡r ❱❛❧✉❡ ▼▲❊ ✭❙❊✮ β ✶✵ ✶ ✶ . ✵✵ ( ✵ . ✶✷ ) ✶ . ✵✵ ( ✵ . ✶✸ ) β ✶✶ − ✵ . ✺ − ✵ . ✹✷ ( ✵ . ✶✽ ) − ✵ . ✹✶ ( ✵ . ✶✽ ) β ✷✵ ✵ . ✺ ✵ . ✹✾ ( ✵ . ✶✸ ) ✵ . ✺✵ ( ✵ . ✶✹ ) β ✷✶ ✵ . ✺ ✵ . ✺✷ ( ✵ . ✶✽ ) ✵ . ✺✸ ( ✵ . ✶✾ ) † ❇❛s❡❞ ♦♥ P♦❧s♦♥ ❡t ❛❧✳ Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r✳ ✺✻ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ P✐❧❧♦✇ ❛♥❞ ❙❝♦tt ✭✷✵✶✷✮ ❡①t❡♥❞ t❤❡ Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r t♦ ♥❡❣❛t✐✈❡ ❜✐♥♦♠✐❛❧ ✭◆❇✮ r❡❣r❡ss✐♦♥ s❡tt✐♥❣ ❈♦♥s✐❞❡r t❤❡ ❢♦❧❧♦✇✐♥❣ ♠♦❞❡❧ ❢♦r ❛ ❝♦✉♥t r✳✈✳ Y i ✿ Γ( y i + r ) d Γ( r ) y i ! ( ✶ − ψ i ) r ψ y i p ( y i | r , β ) = i , r > ✵ , ✇❤❡r❡ � � ❡①♣ ① T i β ❡①♣ ( η i ) � � = ψ i = ① T ✶ + ❡①♣ ( η i ) ✶ + ❡①♣ i β ◆♦t❡ t❤❛t t❤❡ ◆❇ ♣r♦❜❛❜✐❧✐t② ♣❛r❛♠❡t❡r ψ i ✐s ♣❛r❛♠❡t❡r✐③❡❞ ✉s✐♥❣ t❤❡ ❡①♣✐t ❢✉♥❝t✐♦♥ ❚❤✐s ❛❧❧♦✇s ✉s t♦ ❛♣♣❧② t❤❡ s❛♠❡ ♣r♦♣❡rt✐❡s ♦❢ t❤❡ Pó❧②❛✲●❛♠♠❛ ❞❡♥s✐t② ❛s ✐♥ t❤❡ ❧♦❣✐st✐❝ ❝❛s❡ ✺✼ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ ❚❤❡ ♠❡❛♥ ❛♥❞ ✈❛r✐❛♥❝❡ ♦❢ Y i ❛r❡ r ψ i ❊ ( Y i | r , β ) = = r ❡①♣ ( η i ) = µ i ✶ − ψ i r ψ i ❱❛r ( Y i | r , β ) = ( ✶ − ψ i ) ✷ = r ❡①♣ ( η i ) [ ✶ + ❡①♣ ( η i )] = µ i ( ✶ + µ i / r ) ❚❤❡ ♣❛r❛♠❡t❡r α = ✶ / r ❝❛♣t✉r❡s t❤❡ ♦✈❡r❞✐s♣❡rs✐♦♥ ✐♥ t❤❡ ❞❛t❛✱ s✉❝❤ t❤❛t ❛s α → ∞ ✱ t❤❡ ❝♦✉♥ts ❜❡❝♦♠❡ ✐♥❝r❡❛s✐♥❣❧② ❞✐s♣❡rs❡❞ r❡❧❛t✐✈❡ t♦ t❤❡ P♦✐ss♦♥ ❞✐str✐❜✉t✐♦♥ ✺✽ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ ❊①♣❧♦✐t✐♥❣ t❤❡ ❡❛r❧✐❡r ♣r♦♣❡rt② ♦❢ t❤❡ Pó❧②❛✲●❛♠♠❛ ❞✐str✐❜✉t✐♦♥✱ ✐t ❢♦❧❧♦✇s t❤❛t � ∞ ❡ − ω i η ✷ ❡ κ i η i i / ✷ p ( ω i | r + y i , ✵ ) ❞ ω i , p ( Y i | r , β ) ∝ ✵ ✇❤❡r❡ κ i = ( y i − r ) / ✷ ❛♥❞ t❤❡ ω i ✬s ❛r❡ ✐♥❞❡♣❡♥❞❡♥t❧② ❞✐str✐❜✉t❡❞ ❛❝❝♦r❞✐♥❣ t♦ P● ( y i + r , η i ) ✳ ✺✾ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ ❋♦❧❧♦✇✐♥❣ P✐❧❧♦✇ ❛♥❞ ❙❝♦tt✱ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r β ✐s � � − ✶ ✷ ( ③ − ❳ β ) T ❲ ( ③ − ❳ β ) p ( β | ② , r , ω ) ∝ π ( β ) ❡①♣ , ✇❤❡r❡ • ③ ✐s ❛♥ n × ✶ ✈❡❝t♦r ✇✐t❤ ❡❧❡♠❡♥ts z i = y i − r ✷ ω i • ❲ = ❞✐❛❣ ( ω i ) ✻✵ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ ●✐✈❡♥ ❝✉rr❡♥t ✈❛❧✉❡s ❢♦r β ❛♥❞ r ✱ t❤❡ ●✐❜❜s s❛♠♣❧❡r ❢♦r t❤❡ ◆❇ ♠♦❞❡❧ ♣r♦❝❡❡❞s ❛s ❢♦❧❧♦✇s✿ ✶ ❋♦r i = ✶ , . . . , n ✱ ❞r❛✇ ω i ❢r♦♠ ✐ts P● ( y i + r , η i ) ❞✐str✐❜✉t✐♦♥✱ ✇❤❡r❡ η i = ① T i β ✷ ❋♦r i = ✶ , . . . , n ✱ ❞❡✜♥❡ z i = y i − r ✷ ω i � � β ✵ , ❚ − ✶ ✸ ❆ss✉♠✐♥❣ ❛ ◆ p ♣r✐♦r✱ ✉♣❞❛t❡ β ❢r♦♠ ✐ts ◆ p ( ♠ , ❱ ) ❢✉❧❧ ✵ ❝♦♥❞✐t✐♦♥❛❧✱ ✇❤❡r❡ � � − ✶ ❚ ✵ + ❳ T ❲ ❳ ❱ = � � ❚ ✵ β ✵ + ❳ T ❲ ③ ♠ = ❱ ✹ ❯♣❞❛t❡ r ✉s✐♥❣ ❛ r❛♥❞♦♠✲✇❛❧❦ ▼❡tr♦♣♦❧✐s✲❍❛st✐♥❣s st❡♣ ✇✐t❤ ❛ ③❡r♦✲tr✉♥❝❛t❡❞ ♥♦r♠❛❧ ♣r♦♣♦s❛❧ ❞❡♥s✐t②✳ ❆❧t❡r♥❛t✐✈❡❧②✱ ✉♣❞❛t❡ r ✉s✐♥❣ ❛ ❝♦♥❥✉❣❛t❡ ●❛♠♠❛ ❞✐str✐❜✉t✐♦♥ ❛s ❞❡s❝r✐❜❡❞ ✐♥ ❉❛❞❛♥❡❤ ❡t ❛❧✳ ✭✷✵✶✽✮ ❙❡❡ ◆❇❴▼❍✳r ❛♥❞ ◆❇❴●✐❜❜s✳r ❢♦r ❞❡t❛✐❧s ✻✶ ✴ ✶✻✵

❈♦♥❥✉❣❛t❡ ●✐❜❜s ❯♣❞❛t❡ ❢♦r r ✿ ❙t❡♣ ✶ ❩❤♦✉ ❛♥❞ ❈❛r✐♥ ✭✷✵✶✺✮ ❛♥❞ ❉❛❞❛♥❡❤ ❡t ❛❧✳ ✭✷✵✶✽✮ ❞❡s❝r✐❜❡ ❛ t✇♦✲st❡♣ ❝♦♥❥✉❣❛t❡ ●✐❜❜s ✉♣❞❛t❡ ❢♦r r ❚❤❡ ❛♣♣r♦❛❝❤ ✐♥tr♦❞✉❝❡s ❛ s❛♠♣❧❡ ♦❢ ❧❛t❡♥t ❝♦✉♥ts ✱ l i ✱ ✉♥❞❡r❧②✐♥❣ ❡❛❝❤ ♦❜s❡r✈❡❞ ❝♦✉♥t y i ❈♦♥❞✐t✐♦♥❛❧ ♦♥ y i ❛♥❞ r ✱ l i ❤❛s ❛ ❞✐str✐❜✉t✐♦♥ ❞❡✜♥❡❞ ❜② ❛ ❈❤✐♥❡s❡ r❡st❛✉r❛♥t t❛❜❧❡ ✭❈❘❚✮ ❞✐str✐❜✉t✐♦♥✿ � y i l i = u j j = ✶ � � r u j ∼ ❇❡r♥ . r + j − ✶ ❚❤❡ ♥❛♠❡ ✏❈❤✐♥❡s❡ r❡st❛✉r❛♥t t❛❜❧❡✑ ❞❡r✐✈❡s ❢r♦♠ t❤❡ ❢❛❝t t❤❛t u j = ✶ ✐❢ ❛ ♥❡✇ ❝✉st♦♠❡r s✐ts ❛♥❞ ❛♥ ✉♥♦❝❝✉♣✐❡❞ t❛❜❧❡ ✐♥ ❛ ❈❤✐♥❡s❡ r❡st❛✉r❛♥t ✭❛❝❝♦r❞✐♥❣ t♦ ❛ s♦✲❝❛❧❧❡❞ ✏❈❤✐♥❡s❡ r❡st❛✉r❛♥t ♣r♦❝❡ss✑✮✱ ❛♥❞ l i ✐s t❤❡ t♦t❛❧ ♥✉♠❜❡r ♦❢ ♦❝❝✉♣✐❡❞ t❛❜❧❡s ✐♥ t❤❡ r❡st❛✉r❛♥t ❛❢t❡r y i ❝✉st♦♠❡rs ❙♦ ✐♥ ❙t❡♣ ✶✱ ✇❡ ❞r❛✇ l i ( i = ✶ , . . . , n ) ❛❝❝♦r❞✐♥❣ t♦ t❤✐s ❈❘❚ ❞✐str✐❜✉t✐♦♥ ✻✷ ✴ ✶✻✵

❈♦♥❥✉❣❛t❡ ●✐❜❜s ❯♣❞❛t❡ ❢♦r r ✿ ❙t❡♣ ✷ ■♥ ❙t❡♣ ✷✱ t❤❡ ❛✉t❤♦rs ❡①♣❧♦✐t t❤❡ ❢❛❝t t❤❛t t❤❡ ◆❇ ❞✐str✐❜✉t✐♦♥ ❝❛♥ ❜❡ ❞❡r✐✈❡❞ ❢r♦♠ ❛ r❛♥❞♦♠ ❝♦♥✈♦❧✉t✐♦♥ ♦❢ ❧♦❣❛r✐t❤♠✐❝ ❘❱s ❙♣❡❝✐✜❝❛❧❧②✱ t❤❡② ♥♦t❡ t❤❛t✱ ❝♦♥❞✐t✐♦♥❛❧ ♦♥ r ❛♥❞ ψ i ✱ ind l i ∼ P♦✐ [ − r ❧♥ ( ✶ − ψ i )] , ✇❤❡r❡ � � ❡①♣ ① T i β ψ i = � � , i = ✶ , . . . , n . ✶ + ❡①♣ ① T i β ❙❡❡ ❉❛❞❛♥❡❤ ❡t ❛❧✳ ✭✷✵✶✽✮ ❛♥❞ ❩❤♦✉ ❛♥❞ ❈❛r✐♥ ✭✷✵✶✺✮ ❢♦r ❞❡t❛✐❧s ❚❤✉s✱ ✐❢ ✇❡ ❛ss✉♠❡ ❛ ●❛ ( a , b ) ♣r✐♦r ❢♦r r ✱ t❤❡♥ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r r ✐♥ ❙t❡♣ ✷ ✐s � � � n � n r | ❧ , ψ ∼ ●❛ a + l i , b − ❧♥ ( ✶ − ψ i ) , i = ✶ i = ✶ ✇❤❡r❡ ❧ = ( l ✶ , . . . , l n ) T ❛♥❞ ψ = ( ψ ✶ , . . . , ψ n ) T ✳ ❚❤❡ ●✐❜❜s ✉♣❞❛t❡ ✜rst ❞r❛✇s l i ( i = ✶ , . . . , n ) ✐♥❞❡♣❡♥❞❡♥t❧② ❢r♦♠ ❛ ❈❘❚ ❞✐str✐❜✉t✐♦♥✱ ❛♥❞ t❤❡♥ r ❢r♦♠ ✐ts ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ●❛♠♠❛ ❞✐str✐❜✉t✐♦♥ ❣✐✈❡♥ ❧ ❛♥❞ ψ ✻✸ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ ✇✐t❤ ▼❍ ❯♣❞❛t❡ ❢♦r r Hybrid Gibbs-MH Sampler for Negative Binomial Model # Priors (diffuse prior for r) beta0<- rep (0,p) T0<- diag (.01,p) s<-0.01 # Proposal variance -- NOTE: may need to lower this as n_i increases # Inits and Store beta<- rep (0,p) Acc<-0 # MH Acceptance counter ######## # MCMC # ######## for (i in 1:nsim){ # Update r eta<-X%*%beta q<-1/(1+ exp (eta)) # dnbinom fn uses q=1-psi rnew<- rtnorm (1,r, sqrt (s),lower=0) # Proposal ratio<- sum ( dnbinom (y,rnew,q,log=T))- sum ( dnbinom (y,r,q,log=T))+ # Likelihood -- diffuse prior for r dtnorm (rnew,r, sqrt (s),0,log=T)- dtnorm (r,rnew, sqrt (s),0,log=T) # Proposal not symmetric if ( log ( runif (1))<ratio) { r<-rnew Acc<-Acc+1 } # Update beta w<- rpg (n,y+r,eta) # Polya weights z<-(y-r)/(2*w) # Latent response v<- solve ( crossprod (X* sqrt (w))+T0) m<-v%*%(T0%*%beta0+ t ( sqrt (w)*X)%*%( sqrt (w)*z)) beta<- c ( rmvnorm (1,m,v)) # Store if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta R[j]<-r } if (i%%100==0) print (i) # 11 seconds to run 2000 iterations with n=1000 } ✻✹ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r ◆❇ ▼♦❞❡❧ ✇✐t❤ ●✐❜❜s ❯♣❞❛t❡ ❢♦r r Negative Binomial Sampler with Gibbs Update for r beta0<- rep (0,p) T0<- diag (.01,p) a<-b<-0.01 # Gamma hyperparms for r # Inits and Store beta<- rep (0,p) l<- rep (0,n) # Latent counts r<-1 ######## # MCMC # ######## for (i in 1:nsim){ # Update latent counts, l, using CRT distribution for(j in 1:n) l[j]<- sum ( rbinom (y[j],1, round (r/(r+1:y[j]-1),6))) # Could try to avoid loop # Rounding avoids numerical instability # Update r from conjugate gamma distribution given l and beta eta<-X%*%beta psi<- exp (eta)/(1+ exp (eta)) r<- rgamma (1,a+ sum (l),b- sum ( log (1-psi))) # Update beta w<- rpg (n,y+r,eta) # Polya weights z<-(y-r)/(2*w) # Latent response v<- solve ( crossprod (X* sqrt (w))+T0) m<-v%*%(T0%*%beta0+ t ( sqrt (w)*X)%*%( sqrt (w)*z)) beta<- c ( rmvnorm (1,m,v)) # Store if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta R[j]<-r } if (i%%100==0) print (i) # 21 seconds to run 2000 iterations with n=1000 } ✻✺ ✴ ✶✻✵

❊①❛♠♣❧❡ ❚❤❡ ♣r♦❣r❛♠ ◆❇✳r ✜ts t❤❡ ❢♦❧❧♦✇✐♥❣ ◆❇ ♠♦❞❡❧✿ Γ( y i + r ) d Γ( r ) y i ! ( ✶ − ψ i ) r ψ y i p ( y i | r , β ) = i , r > ✵ , ✇❤❡r❡ ❡①♣ ( β ✵ + β ✶ x i ) ψ i = ✶ + ❡①♣ ( β ✵ + β ✶ x i ) ❚❤❡ r❡s✉❧ts✱ ❜❛s❡❞ ♦♥ ✷✵✵✵ ✐t❡r❛t✐♦♥s ✇✐t❤ ❛ ❜✉r♥✲✐♥ ♦❢ ✶✵✵✵✱ ❛r❡✿ ❚❛❜❧❡ ✼✿ ❘❡s✉❧ts ❢♦r ◆❇ ▼♦❞❡❧✳ ❚r✉❡ P♦st❡r✐♦r P♦st❡r✐♦r ▼❡❛♥ ✭❙❉✮ † ▼❡❛♥ ✭❙❉✮ ‡ P❛r❛♠❡t❡r ❱❛❧✉❡ ▼▲❊ ✭❙❊✮ β ✶ ✶ . ✵ ✵ . ✾✸ ( ✵ . ✵✽ ) ✵ . ✾✷ ( ✵ . ✵✽ ) ✵ . ✾✷ ( ✵ . ✵✼ ) β ✷ ✵ . ✺ ✵ . ✺✵ ( ✵ . ✵✹ ) ✵ . ✺✵ ( ✵ . ✵✹ ) ✵ . ✺✵ ( ✵ . ✵✹ ) ✶ . ✵ ✶ . ✶✵ ( ✵ . ✵✽ ) ✶ . ✶✵ ( ✵ . ✵✼ ) ✶ . ✶✶ ( ✵ . ✵✼ ) r † Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ▼❍ ✉♣❞❛t❡ ❢♦r r ✳ ❆❝❝❡♣t❛♥❝❡ r❛t❡ ❂ ✸✼✪✳ ‡ Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ●✐❜❜s ✉♣❞❛t❡ ❢♦r r ✳ ✻✻ ✴ ✶✻✵

❚r❛❝❡ P❧♦ts ❢♦r ◆❇ ▼♦❞❡❧ ✇✐t❤ ▼❍ ❯♣❞❛t❡ ❢♦r r 2.2 2.0 β 1 1.8 1.6 0 200 400 600 800 1000 Iteration 0.8 0.6 β 2 0.4 0.2 0 200 400 600 800 1000 Iteration 1.4 1.2 r 1.0 0.8 0 200 400 600 800 1000 Iteration ✻✼ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ❩❡r♦✲■♥✢❛t❡❞ ❈♦✉♥t ❉❛t❛ ❩❡r♦✲✐♥✢❛t❡❞ ❝♦✉♥t ❞❛t❛ ❛r✐s❡ ✇❤❡♥ t❤❡ ❞❛t❛ ❝♦♥t❛✐♥ ❛ ❧❛r❣❡r ♣r♦♣♦rt✐♦♥ ♦❢ ③❡r♦s t❤❛♥ ♣r❡❞✐❝t❡❞ ❜② ❛♥ ♦r❞✐♥❛r② ❝♦✉♥t ♠♦❞❡❧ s✉❝❤ ❛s t❤❡ ◆❇ ❩❡r♦✲✐♥✢❛t❡❞ ♠♦❞❡❧s ❤❛✈❡ ❜❡❡♥ ♣r♦♣♦s❡❞ t♦ ❛❞❞r❡ss t❤❡ ♦✈❡r✲❛❜✉♥❞❛♥❝❡ ♦❢ ③❡r♦s ❩❡r♦✲✐♥✢❛t❡❞ ♠♦❞❡❧s ❛r❡ ♠✐①t✉r❡s ♦❢ ❛ ♣♦✐♥t ♠❛ss ❛t ③❡r♦✱ r❡♣r❡s❡♥t✐♥❣ t❤❡ ❡①❝❡ss ③❡r♦s✱ ❛♥❞ ❛ ❝♦✉♥t ❞✐str✐❜✉t✐♦♥ ❢♦r t❤❡ r❡♠❛✐♥✐♥❣ ✈❛❧✉❡s ❇② ❝♦♥str✉❝t✐♦♥✱ ③❡r♦✲✐♥✢❛t❡❞ ♠♦❞❡❧s ♣❛rt✐t✐♦♥ ③❡r♦s ✐♥t♦ t✇♦ t②♣❡s✿ • ❙tr✉❝t✉r❛❧ ③❡r♦s ❝♦rr❡s♣♦♥❞✐♥❣ t♦ ✐♥❞✐✈✐❞✉❛❧s ✇❤♦ ❛r❡ ♥♦t ❛t r✐s❦ ❢♦r ❛♥ ❡✈❡♥t✱ ❛♥❞ t❤❡r❡❢♦r❡ ❤❛✈❡ ♥♦ ♦♣♣♦rt✉♥✐t② ❢♦r ❛ ♣♦s✐t✐✈❡ ❝♦✉♥t • ❆t✲r✐s❦ ③❡r♦s✱ ✇❤✐❝❤ ❛♣♣❧② t♦ ❛ ❧❛t❡♥t ❝❧❛ss ♦❢ ✐♥❞✐✈✐❞✉❛❧s ✇❤♦ ❛r❡ ❛t r✐s❦ ❢♦r ❛♥ ❡✈❡♥t ❜✉t ♥❡✈❡rt❤❡❧❡ss ❤❛✈❡ ❛♥ ♦❜s❡r✈❡❞ r❡s♣♦♥s❡ ♦❢ ③❡r♦ ✻✽ ✴ ✶✻✵

❩❡r♦✲■♥✢❛t❡❞ ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ ❈♦♥s✐❞❡r t❤❡ ❢♦❧❧♦✇✐♥❣ ③❡r♦✲✐♥✢❛t❡❞ ♥❡❣❛t✐✈❡ ❜✐♥♦♠✐❛❧ ✭❩■◆❇✮ ♠♦❞❡❧ ( ✶ − φ i ) + φ i ( ✶ − ψ i ) r Pr ( Y i = ✵ ) = Γ( y i + r ) Γ( r ) y i ! ( ✶ − ψ i ) r ψ y i Pr ( Y i = y i ) = φ i i , y i = ✶ , ✷ , . . . � � ① T ❡①♣ i β ❡①♣ ( η i ) � � = ψ i = ✶ + ❡①♣ ① T ✶ + ❡①♣ ( η i ) i β ❚❤❡ ♣❛r❛♠❡t❡r φ i ❞❡♥♦t❡s t❤❡ ♣r♦❜❛❜✐❧✐t② ♦❢ ❜❡❧♦♥❣✐♥❣ t♦ t❤❡ ✏❛t✲r✐s❦✑ ❝❧❛ss ✶ − φ i ❞❡♥♦t❡s t❤❡ ♣r♦❜❛❜✐❧✐t② ♦❢ ❛♥ ❡①❝❡ss ✭✐✳❡✳✱ str✉❝t✉r❛❧✮ ③❡r♦ ✻✾ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ❢♦r ❩❡r♦✲■♥✢❛t❡❞ ❈♦✉♥t ❉❛t❛ ❲❡ ❝❛♥ r❡✇r✐t❡ t❤❡ ❩■◆❇ ♠♦❞❡❧ ❜② ✐♥tr♦❞✉❝✐♥❣ ❛ ❧❛t❡♥t ✏❛t✲r✐s❦✑ ✐♥❞✐❝❛t♦r ✈❛r✐❛❜❧❡✱ W i ✱ Y i ∼ ( ✶ − φ i ) ✶ ( W i = ✵ ∧ Y i = ✵ ) + φ i ◆❇ ( µ i , r ) ✶ ( W i = ✶ ) , i = ✶ , . . . , n ❈♦♠♠❡♥ts✿ • ❲✐t❤ ♣r♦❜❛❜✐❧✐t② ✶ − φ i ✱ W i = Y i = ✵ ✭str✉❝t✉r❛❧ ③❡r♦✮ • ❲✐t❤ ♣r♦❜❛❜✐❧✐t② φ i ✱ W i = ✶ ❛♥❞ Y i ✐s ❞r❛✇♥ ❢r♦♠ ❛ ◆❇ ❞✐str✐❜✉t✐♦♥ ✇✐t❤ ♠❡❛♥ µ i ❛♥❞ ❞✐s♣❡rs✐♦♥ ♣❛r❛♠❡t❡r r > ✵ • µ i = ❊ ( Y i | β , r , W i = ✶ ) ✐s t❤❡ ♠❡❛♥ ❝♦✉♥t ❛♠♦♥❣ t❤♦s❡ ✐♥ t❤❡ ❛t✲r✐s❦ ❝❧❛ss ✭❝♦♥❞✐t✐♦♥❛❧ ♦♥ W i ❂✶✮ • r ❝❛♣t✉r❡s ♦✈❡r❞✐s♣❡rs✐♦♥ ✐♥ t❤❡ ❛t✲r✐s❦ ❝❧❛ss • ❚❤❡ ♦✈❡r❛❧❧ ✭✉♥❝♦♥❞✐t✐♦♥❛❧✮ ♠❡❛♥ ✐s ❊ ( Y i | β , r ) = φ i µ i • ■❢ φ i = ✶✱ t❤❡ ❩■◆❇ r❡❞✉❝❡s t♦ ❛♥ ◆❇ ♠♦❞❡❧ • ✐❢ φ i = ✵✱ ✇❡ ❤❛✈❡ ❛ ♣♦✐♥t ♠❛ss ❛t ③❡r♦ ✼✵ ✴ ✶✻✵

❩❡r♦✲■♥✢❛t❡❞ ◆❡❣❛t✐✈❡ ❇✐♥♦♠✐❛❧ ▼♦❞❡❧ ❚②♣✐❝❛❧❧②✱ ✇❡ ♠♦❞❡❧ φ i ✉s✐♥❣ ❛ ❧♦❣✐t ♠♦❞❡❧ ❛♥❞ Y i | W i = ✶ ✉s✐♥❣ ❛ ◆❇ ♠♦❞❡❧✿ ❧♦❣✐t [ Pr ( W i = ✶ | α )] = ① T ❧♦❣✐t ( φ i ) = i α = η ✶ i Γ( y i + r ) d Γ( r ) y i ! ( ✶ − ψ i ) r ψ y i p ( y i | r , β , W i = ✶ ) = ∀ i s✳t✳ W i = ✶ , ✇❤❡r❡ i � � ① T ❡①♣ i β ❡①♣ ( η ✷ i ) � � = ψ i = ✶ + ❡①♣ ( η ✷ i ) . ① T ✶ + ❡①♣ i β ✼✶ ✴ ✶✻✵

●✐❜❜s ❙❛♠♣❧❡r ❢♦r ❩■◆❇ ▼♦❞❡❧ ❚❤✐s s✉❣❣❡sts t❤❡ ❢♦❧❧♦✇✐♥❣ ●✐❜❜s s❛♠♣❧❡r✿ ✶ ●✐✈❡♥ ❝✉rr❡♥t ♣❛r❛♠❡t❡r ✈❛❧✉❡s✱ ✉♣❞❛t❡ α ✉s✐♥❣ t❤❡ ●✐❜❜s s❛♠♣❧❡r ♣r♦♣♦s❡❞ ❜② P♦❧s♦♥ ❡t ❛❧✳ ❢♦r ❧♦❣✐st✐❝ r❡❣r❡ss✐♦♥ ✷ ❈♦♥❞✐t✐♦♥❛❧ ♦♥ W i = ✶✱ ✉♣❞❛t❡ β ✉s✐♥❣ t❤❡ ◆❇ ●✐❜❜s s❛♠♣❧❡r ♣r♦♣♦s❡❞ ❜② P✐❧❧♦✇ ❛♥❞ ❙❝♦tt ✸ ❯♣❞❛t❡ r ✉s✐♥❣ ❛ r❛♥❞♦♠✲✇❛❧❦ ▼❡tr♦♣♦❧✐s✲❍❛st✐♥❣s st❡♣ ♦r ✉s✐♥❣ ❛ ❝♦♥❥✉❣❛t❡ ●❛♠♠❛ ✉♣❞❛t❡ ❛s ✐♥ ❉❛❞❡♥❡❤ ❡t ❛❧✳ ✭✷✵✶✽✮ ✹ ❯♣❞❛t❡ t❤❡ ❧❛t❡♥t ❛t✲r✐s❦ ✐♥❞✐❝❛t♦rs✱ W ✶ , . . . , W n ✱ ❢r♦♠ t❤❡✐r ❞✐s❝r❡t❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❞✐str✐❜✉t✐♦♥s ❚❤❡ ♦♥❧② ♥❡✇ st❡♣ ✐s st❡♣ ✭✹✮✱ t❤❡ ✉♣❞❛t❡ ❢♦r W i ✼✷ ✴ ✶✻✵

❯♣❞❛t❡ ❢♦r W i ❚❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r W i ✐s ❛ ❞✐s❝r❡t❡ ❞✐str✐❜✉t✐♦♥ ✇✐t❤ ♣r♦❜❛❜✐❧✐t✐❡s t❤❛t ❞❡♣❡♥❞ ♦♥ ✇❤❡t❤❡r t❤❡ ♦❜s❡r✈❡❞ ❝♦✉♥t✱ y i ✱ ✐s ③❡r♦ ♦r ♥♦♥✲③❡r♦ ■❢ y i > ✵✱ t❤❡♥ s✉❜❥❡❝t i ❜❡❧♦♥❣s t♦ t❤❡ ❛t✲r✐s❦ ❝❧❛ss✱ ❛♥❞ ❤❡♥❝❡ ❜② ❞❡✜♥✐t✐♦♥✱ W i = ✶ ✇✐t❤ ♣r♦❜❛❜✐❧✐t② ✶ ❈♦♥✈❡rs❡❧②✱ ✐❢ y i = ✵✱ t❤❡♥ ✇❡ ♦❜s❡r✈❡ ❡✐t❤❡r ❛ str✉❝t✉r❛❧ ③❡r♦ ✭✐♠♣❧②✐♥❣ t❤❛t W i = ✵✮ ♦r ❛♥ ❛t✲r✐s❦ ③❡r♦ ✭✐♠♣❧②✐♥❣ W i = ✶✮ ❙❡❡ ◆❡❡❧♦♥ ✭✷✵✶✽✮ ❢♦r ❞❡t❛✐❧s ✼✸ ✴ ✶✻✵

❯♣❞❛t❡ ❢♦r W i ∗ ■♥ t❤❡ ❧❛tt❡r ❝❛s❡✱ ✇❡ ❞r❛✇ W i ❢r♦♠ ❛ ❇❡r♥♦✉❧❧✐ ❞✐str✐❜✉t✐♦♥ ✇✐t❤ ♣r♦❜❛❜✐❧✐t② θ i = Pr ( W i = ✶ | y i = ✵ , r❡st ) = Pr ( ❛t✲r✐s❦ ③❡r♦ | ❛t r✐s❦ ♦r str✉❝t✉r❛❧ ③❡r♦ ) φ i υ r i = i ) , ✶ − φ i ( ✶ − υ r ✇❤❡r❡ • φ i = ❡①♣ ( η ✶ i ) / [ ✶ + ❡①♣ ( η ✶ i )] ✐s t❤❡ ✉♥❝♦♥❞✐t✐♦♥❛❧ ♣r♦❜❛❜✐❧✐t② t❤❛t W i = ✶ • υ i = ✶ − ψ i ✱ ✇❤❡r❡ ψ i ✐s t❤❡ ◆❇ ❡✈❡♥t ♣r♦❜❛❜✐❧✐t② ❙❡❡ ❩■◆❇✳r ❢♦r ❞❡t❛✐❧s ✼✹ ✴ ✶✻✵

❊①❛♠♣❧❡ ❉❛t❛ ❍✐st♦❣r❛♠ 60 AIC for ZINB is 104 points lower than for NB model 50 40 Percent 30 20 10 0 0 2 4 6 8 11 14 17 20 23 26 29 32 35 38 44 60 Count ✼✺ ✴ ✶✻✵

❍②❜r✐❞ ●✐❜❜s✲▼❍ ❆❧❣♦r✐t❤♠ ❢♦r ❩■◆❇ ▼♦❞❡❧ Hybrid Gibbs-MH Sampler for the ZINB Model for (i in 1:nsim){ # Update alpha mu<-X%*%alpha w<- rpg (n,1,mu) z<-(y1-1/2)/w # Latent response v<- solve ( crossprod (X* sqrt (w))+T0a) m<-v%*%(T0a%*%alpha0+ t (w*X)%*%z) alpha<- c ( rmvnorm (1,m,v)) # Update latent class indicator y1 (= W in slides) eta<-X%*%alpha phi<- exp (eta)/(1+ exp (eta)) # At-risk probability theta<-phi*(q^r)/(phi*(q^r)+1-phi) # Cond prob that y1=1 given y=0 -- i.e. Pr(chance zero|observed y1[y==0]<- rbinom (n0,1,theta[y==0]) # If y=0, draw "chance zero" w.p. theta; if y=1, then y1=1 n1<- sum (y1) # Update beta conditional on y1=1 eta1<-X[y1==1,]%*%beta w<- rpg (n1,y[y1==1]+r,eta1) # Polya weights z<-(y[y1==1]-r)/(2*w) # Latent response v<- solve ( crossprod (X[y1==1,]* sqrt (w))+T0b) m<-v%*%(T0b%*%beta0+ t ( sqrt (w)*X[y1==1,])%*%( sqrt (w)*z)) beta<- c ( rmvnorm (1,m,v)) eta<-X%*%beta q<-1/(1+ exp (eta)) # Update r rnew<- rtnorm (1,r, sqrt (s),lower=0) ratio<- sum ( dnbinom (y[y1==1],rnew,q[y1==1],log=T))- sum ( dnbinom (y[y1==1],r,q[y1==1],log=T))+ # Diffuse prior dtnorm (r,rnew, sqrt (s),0,log=T) - dtnorm (rnew,r, sqrt (s),0,log=T) if ( log ( runif (1))<ratio) { r<-rnew Acc<-Acc+1 } # Store if (i> burn & i%%thin==0) { j<-(i-burn)/thin Alpha[j,]<-alpha Beta[j,]<-beta R[j]<-r } if (i%%100==0) print (i) # 11 seconds to run 2000 iterations with n=1000 } ✼✻ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r ❩■◆❇ ▼♦❞❡❧ ✇✐t❤ ●✐❜❜s ❯♣❞❛t❡ ❢♦r r Gibbs Sampler for ZINB Model for (i in 1:nsim){ # Update alpha mu<-X%*%alpha w<- rpg (n,1,mu) z<-(y1-1/2)/w # Latent response v<- solve ( crossprod (X* sqrt (w))+T0a) m<-v%*%(T0a%*%alpha0+ t (w*X)%*%z) alpha<- c ( rmvnorm (1,m,v)) # Update latent class indicator y1 (= W in slides) eta<-X%*%alpha phi<- exp (eta)/(1+ exp (eta)) # At-risk probability theta<-phi*(q^r)/(phi*(q^r)+1-phi) # Cond prob that y1=1 given y=0 -- i.e. Pr(chance zero|observed y1[y==0]<- rbinom (n0,1,theta[y==0]) # If y=0, draw "chance zero" w.p. theta; if y=1, then y1=1 n1<- sum (y1) # Update beta conditional on y1=1 eta1<-X[y1==1,]%*%beta w<- rpg (n1,y[y1==1]+r,eta1) # Polya weights z<-(y[y1==1]-r)/(2*w) # Latent response v<- solve ( crossprod (X[y1==1,]* sqrt (w))+T0b) m<-v%*%(T0b%*%beta0+ t ( sqrt (w)*X[y1==1,])%*%( sqrt (w)*z)) beta<- c ( rmvnorm (1,m,v)) eta<-X%*%beta q<-1/(1+ exp (eta)) # Update latent counts, l l<- rep (0,n1) ytmp<-y[y1==1] for(j in 1:n1) l[j]<- sum ( rbinom (ytmp[j],1, round (r/(r+1:ytmp[j]-1),6))) # Update r from conjugate gamma distribution given l and psi eta<-X[y1==1,]%*%beta psi<- exp (eta)/(1+ exp (eta)) r<- rgamma (1,a+ sum (l),b- sum ( log (1-psi))) # Store if (i> burn & i%%thin==0) { j<-(i-burn)/thin Alpha[j,]<-alpha Beta[j,]<-beta R[j]<-r } if (i%%100==0) print (i) # 15 seconds to run 2000 iterations with n=1000 } ✼✼ ✴ ✶✻✵

❊①❛♠♣❧❡ ❚❤❡ ♣r♦❣r❛♠ ❩■◆❇✳r ✜ts t❤❡ ❢♦❧❧♦✇✐♥❣ ❩■◆❇ ♠♦❞❡❧✿ ❧♦❣✐t ( φ i ) = α ✶ + α ✷ x i Γ( y i + r ) d Γ( r ) y i ! ( ✶ − ψ i ) r ψ y i p ( y i | r , β , W i = ✶ ) = ∀ i s✳t✳ W i = ✶ , ✇❤❡r❡ i ❡①♣ ( β ✶ + β ✷ x i ) ψ i = ✶ + ❡①♣ ( β ✶ + β ✷ x i ) ❚❤❡ r❡s✉❧ts✱ ❜❛s❡❞ ♦♥ ✷✵✵✵ ✐t❡r❛t✐♦♥s ✇✐t❤ ❛ ❜✉r♥✲✐♥ ♦❢ ✶✵✵✵✱ ❛r❡✿ ❚❛❜❧❡ ✽✿ ❘❡s✉❧ts ❢♦r ❩■◆❇ ▼♦❞❡❧✳ ❚r✉❡ P♦st❡r✐♦r P♦st❡r✐♦r ▼❡❛♥ ✭❙❉✮ † ▼❡❛♥ ✭❙❉✮ ‡ P❛r❛♠❡t❡r ❱❛❧✉❡ ▼▲❊ ✭❙❊✮ α ✶ − ✵ . ✺ − ✵ . ✻✸ ( ✵ . ✶✶ ) − ✵ . ✻✷ ( ✵ . ✶✵ ) − ✵ . ✻✹ ( ✵ . ✶✶ ) α ✷ ✵ . ✺ ✵ . ✺✻ ( ✵ . ✶✹ ) ✵ . ✺✺ ( ✵ . ✶✹ ) ✵ . ✺✼ ( ✵ . ✶✹ ) β ✶ ✷ . ✵ ✶ . ✾✹ ( ✵ . ✶✸ ) ✶ . ✾✻ ( ✵ . ✶✷ ) ✶ . ✾✹ ( ✵ . ✶✸ ) β ✷ ✵ . ✺ ✵ . ✸✽ ( ✵ . ✶✶ ) ✵ . ✸✾ ( ✵ . ✶✶ ) ✵ . ✸✾ ( ✵ . ✶✶ ) ✶ . ✵ ✶ . ✶✶ ( ✵ . ✶✸ ) ✶ . ✶✵ ( ✵ . ✶✶ ) ✶ . ✶✸ ( ✵ . ✶✸ ) r † Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ▼❍ ✉♣❞❛t❡ ❢♦r r ✳ ❆❝❝❡♣t❛♥❝❡ r❛t❡ ❂ ✹✵✪✳ ‡ Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ●✐❜❜s ✉♣❞❛t❡ ❢♦r r ✳ ✼✽ ✴ ✶✻✵

❚r❛❝❡ P❧♦ts ❢♦r ▼❍ ❩■◆❇ ▼♦❞❡❧ −0.3 1.0 −0.5 0.6 α 1 α 2 −0.7 0.2 −0.9 0 200 400 600 800 1000 0 200 400 600 800 1000 1:lastit 1:lastit 0.8 2.2 0.6 2.0 β 1 β 2 0.4 1.8 0.2 1.6 0 200 400 600 800 1000 0 200 400 600 800 1000 1:lastit 1:lastit 1.4 1.2 r 1.0 0.8 0 200 400 600 800 1000 1:lastit ✼✾ ✴ ✶✻✵

❚r❛❝❡ P❧♦ts ❢♦r ●✐❜❜s ❩■◆❇ ▼♦❞❡❧ −0.2 0.8 −0.6 α 1 α 2 0.4 −1.0 0.0 0 200 400 600 800 1000 0 200 400 600 800 1000 1:lastit 1:lastit 0.8 2.2 0.6 2.0 β 1 β 2 0.4 1.8 0.2 1.6 0 200 400 600 800 1000 0 200 400 600 800 1000 1:lastit 1:lastit 1.4 1.2 r 1.0 0.8 0 200 400 600 800 1000 1:lastit ✽✵ ✴ ✶✻✵

❇❛②❡s✐❛♥ ▼♦❞❡❧✐♥❣ ❙tr❛t❡❣✐❡s ❢♦r ●❡♥❡r❛❧✐③❡❞ ▲✐♥❡❛r ▼♦❞❡❧s✱ P❛rt ✷ ❘❡❛❞✐♥❣ ✿ ❍♦✛ ❙❡❝t✐♦♥ ✼✳✸❀ ❲❛❦❡✜❡❧❞ ❙❡❝t✐♦♥s ✶✵✳✺✳✶✱✶✶✳✷✳✽✱ ✶✶✳✷✳✾ ✭P❡♥❛❧✐③❡❞ ❘❡❣r❡ss✐♦♥✮❀ ❲❛❦❡✜❡❧❞ ❙❡❝t✐♦♥s ✽✳✻✱ ✽✳✼ ✭▲✐♥❡❛r ▼✐①❡❞ ▼♦❞❡❧s✮❀ ◆❡❡❧♦♥ ✭✷✵✶✽✮ ❙❡❝t✐♦♥s ✷✳✸✕✸✳✷❀ ❋rü❤✇✐rt❤✲❙❝❤♥❛tt❡r ❛♥❞ P②♥❡ ✭✷✵✶✵✮❀ ◆❡❡❧♦♥ ❡t ❛❧✳ ✭✷✵✶✺✮ ❋❛❧❧ ✷✵✶✽ ✽✶ ✴ ✶✻✵

❇❛②❡s✐❛♥ P❡♥❛❧✐③❡❞ ▲✐♥❡❛r ❙♣❧✐♥❡s ▲❡t✬s ❝♦♥s✐❞❡r t❤❡ ♣r♦❜❧❡♠ ♦❢ ♠♦❞❡❧✐♥❣ ❛ ♠❡❛♥ r❡s♣♦♥s❡ ✉s✐♥❣ ❛ ♣✐❡❝❡✇✐s❡ ❧✐♥❡❛r ✭P❲▲✮ s♣❧✐♥❡ ❢✉♥❝t✐♦♥ ❆s ❛♥ ✐❧❧✉str❛t✐♦♥✱ ❝♦♥s✐❞❡r t❤❡ ❢♦❧❧♦✇✐♥❣ s❝❛tt❡r♣❧♦t s❤♦✇✐♥❣ t❤❡ r❡❧❛t✐♦♥s❤✐♣ ❜❡t✇❡❡♥ ❛❣❡ ❛♥❞ ❧♦❣ ❈✲♣❡♣t✐❞❡ ❝♦♥❝❡♥tr❛t✐♦♥ ❛♠♦♥❣ ✹✸ ❞✐❛❜❡t✐❝ ❝❤✐❧❞r❡♥ ✷ ❚❤✐s ❡①❛♠♣❧❡ ✐s ❢♦r ✐♥❞❡♣❡♥❞❡♥t ❞❛t❛✱ ❜✉t t❤❡ ❛♣♣r♦❛❝❤ ❡❛s✐❧② ❡①t❡♥❞s t♦ r❡♣❡❛t❡❞ ♠❡❛s✉r❡s ❞❛t❛ ❋✐❣✉r❡ ✷✿ ❙❝❛tt❡r♣❧♦t ♦❢ ❧♦❣ ❈✲♣❡♣t✐❞❡ ❝♦♥❝❡♥tr❛t✐♦♥ ❜② ❛❣❡✳ ● 0.8 ● ● ● ● ● ● ● ● ● 0.7 ● ● ● ● ● ● ● Log Concentration ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.6 ● ● ● ● ● ● ● 0.5 ● 0 2 4 6 8 10 12 14 16 Age (Y ears) ✷ ❝✳❢✳✱ ❋✐t③♠❛✉r✐❝❡ ❡t ❛❧✳✱ ❈❤❛♣t❡r ✶✾ ✽✷ ✴ ✶✻✵

❇❛②❡s✐❛♥ P❡♥❛❧✐③❡❞ ▲✐♥❡❛r ❙♣❧✐♥❡s ❚❤❡ ❣♦❛❧ ✐s t♦ ✢❡①✐❜❧② ♠♦❞❡❧ ❧♦❣ ♣❡♣t✐❞❡ ❝♦♥❝❡♥tr❛t✐♦♥ ❜② ❛❣❡ ❆ r❡❧❛t✐✈❡❧② s✐♠♣❧❡ ❝❤♦✐❝❡ ✐s t♦ ✜t t❤❡ ❢♦❧❧♦✇✐♥❣ P❲▲ ♠♦❞❡❧✿ � K ❊ ( Y ij | x ij ) = β ✵ + β ✶ x ij + b k ( x ij − c k ) + , i = ✶ , . . . , n k = ✶ ✇❤❡r❡ c ✶ , . . . , c K ❛r❡ ✐♥t❡r✐♦r ❦♥♦t ❧♦❝❛t✐♦♥s✱ b ✶ , . . . , b K ❛r❡ s♣❧✐♥❡ ❝♦❡✣❝✐❡♥ts✱ ❛♥❞ ❢♦r ❛ ♥✉♠❜❡r u ✱ u + = u ✐❢ u > ✵ ❛♥❞ ✵ ♦t❤❡r✇✐s❡✳ ◆♦t❡ t❤❛t ✐❢ ❛❧❧ b k = ✵✱ ✇❡ r❡❞✉❝❡ t♦ ❛ str❛✐❣❤t✲❧✐♥❡ r❡❣r❡ss✐♦♥ ❢✉♥❝t✐♦♥ ✽✸ ✴ ✶✻✵

❇❛②❡s✐❛♥ P❡♥❛❧✐③❡❞ ▲✐♥❡❛r ❙♣❧✐♥❡s ❋✐tt✐♥❣ t❤❡ ♠♦❞❡❧ ✐♥ ❙❆❙ ✇✐t❤ ✶✵ ❦♥♦ts ❢r♦♠ ❛❣❡s ✺✲✶✹✱ ✇❡ s❡❡ t❤❛t t❤❡ ♠♦❞❡❧ ✐s r❡❛s♦♥❛❜❧② ✢❡①✐❜❧❡ ❜✉t ♥♦t s♠♦♦t❤ ❋✐❣✉r❡ ✸✿ ❊st✐♠❛t❡❞ ♠❡❛♥ r❡❣r❡ss✐♦♥ ❢✉♥❝t✐♦♥ ❢♦r ♣✐❡❝❡✇✐s❡ ❧✐♥❡❛r ♠♦❞❡❧✳ ● 0.8 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.7 ● Log Concentration ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.6 ● ● ● ● ● ● ● 0.5 ● 0 2 4 6 8 10 12 14 16 Age (Y ears) ✽✹ ✴ ✶✻✵

P❡♥❛❧✐③❡❞ ▲✐♥❡❛r ❙♣❧✐♥❡s ❚♦ ✐♠♣♦s❡ s♠♦♦t❤♥❡ss✱ ✇❡ ❝❛♥ ✏s❤r✐♥❦✑ t❤❡ s♣❧✐♥❡ ❝♦❡✣❝✐❡♥ts b ✶ , . . . , b K ❜② ♣❡♥❛❧✐③✐♥❣ ❧❛r❣❡ ✈❛❧✉❡s ❊ss❡♥t✐❛❧❧② ✇❡ s❤r✐♥❦ ❧❛r❣❡ s♣❧✐♥❡ ❝♦❡✣❝✐❡♥ts t♦✇❛r❞ ③❡r♦✱ s♦ t❤❛t t❤❡r❡ ❛r❡ ♥♦ ❡①tr❡♠❡ ✈❛❧✉❡s ❚❤✐s ❝❛♥ ❜❡ ✉s❡❢✉❧ ✐❢ ✇❡ ❤❛✈❡ ❛ ❧❛r❣❡ ♥✉♠❜❡r ♦❢ ❦♥♦ts✱ s✐♥❝❡ t❤✐s ❝❛♥ ✐♥tr♦❞✉❝❡ ❝♦❧❧✐♥❡❛r✐t② ❛♠♦♥❣ t❤❡ P❲▲ ❜❛s✐s ❢✉♥❝t✐♦♥s ❉♦✐♥❣ s♦ ✇♦✉❧❞ ❧❡❛❞ t♦ ❧❛r❣❡ ❙❊s ❢♦r t❤❡ s♣❧✐♥❡ ❝♦❡✣❝✐❡♥ts P❡♥❛❧✐③✐♥❣ t❤❡ ❝♦❡✣❝✐❡♥ts st❛❜✐❧✐③❡s t❤❡ ❡st✐♠❛t❡s ❛♥❞ ❛✈♦✐❞s ♦✈❡r✜tt✐♥❣✱ ✇❤✐❝❤ ❝❛♥ ❛r✐s❡ ✇❤❡♥ ✇❡ ✜t ❛♥ ♦✈❡r❧② ✏♥♦✐s②✑ ❝✉r✈❡ ✽✺ ✴ ✶✻✵

P❡♥❛❧✐③❡❞ ❘❡❣r❡ss✐♦♥ ❘❛t❤❡r t❤❛♥ ♠✐♥✐♠✐③✐♥❣ t❤❡ r❡s✐❞✉❛❧ s✉♠ ♦❢ sq✉❛r❡s � � �� ✷ � n � K Q = Y i − β ✵ + β ✶ x i + b k ( x i − c k ) + , i = ✶ k = ✶ ♣❡♥❛❧✐③❡❞ r❡❣r❡ss✐♦♥ ♠✐♥✐♠✐③❡s Q s✉❜❥❡❝t t♦ ❛ ❝♦♥str❛✐♥t t❤❛t r❡str✐❝ts t❤❡ s✐③❡ ♦❢ t❤❡ { b k } ✱ t❤✉s s♠♦♦t❤✐♥❣ ♦r r❡❣✉❧❛r✐③✐♥❣ t❤❡ s♣❧✐♥❡ ❝✉r✈❡ P♦♣✉❧❛r ❝♦♥str❛✐♥ts ✐♥❝❧✉❞❡ ✸ ✶ ♠❛① | b k | < t ✷ � K k = ✶ | b k | < t ✭▲❛ss♦ ❝♦♥str❛✐♥t✮ ✸ � K k < t = ❜ T ❜ < t ✭❘✐❞❣❡ ❝♦♥str❛✐♥t✮ k = ✶ b ✷ ✹ ❜ T ❉❜ < t ❢♦r s♦♠❡ K × K ♣♦s✐t✐✈❡ s❡♠✐✲❞❡✜♥✐t❡ ♣❡♥❛❧t② ♠❛tr✐① ✱ ❉ ✭❡✳❣✳✱ ❛ s♣❛t✐❛❧ s♠♦♦t❤✐♥❣ ♠❛tr✐①✮ ✸ ❘✉♣♣❡rt ❡t ❛❧✳✱ ♣❛❣❡ ✻✺✳ ✽✻ ✴ ✶✻✵

P❡♥❛❧✐③❡❞ ❘❡❣r❡ss✐♦♥ P❡♥❛❧✐③❡❞ r❡❣r❡ss✐♦♥ ♠✐♥✐♠✐③❡s t❤❡ ♣❡♥❛❧✐③❡❞ s✉♠ ♦❢ sq✉❛r❡s � � �� ✷ � n � K + λ ❜ T ❉❜ , Q ( λ ) = Y i − β ✵ + β ✶ x i + b k ( x i − c k ) + i = ✶ k = ✶ ✇❤❡r❡ • λ ❜ T ❉❜ ✐s ❦♥♦✇♥ ❛s t❤❡ r♦✉❣❤♥❡ss ♣❡♥❛❧t② ♣❡♥❛❧✐③✐♥❣ ♦✈❡r❧② ✏r♦✉❣❤✑ ♦r ♥♦✐s② ❢❡❛t✉r❡s ♦❢ t❤❡ ❝✉r✈❡✱ ❛♥❞ • λ ✐s ❛ t✉♥✐♥❣ ♣❛r❛♠❡t❡r t❤❛t ❝♦♥tr♦❧s t❤❡ ❞❡❣r❡❡ ♦❢ s♠♦♦t❤♥❡ss✱ ✇✐t❤ ✐♥❝r❡❛s✐♥❣ s♠♦♦t❤♥❡ss ❛s λ → ∞ λ ❝❛♥ ❜❡ ✈✐❡✇❡❞ ❛s ❛ ▲❛❣r❛♥❣❡ ♠✉❧t✐♣❧✐❡r ✱ ✇❤✐❝❤ ✐s ❛ t❡❝❤♥✐q✉❡ ❝♦♠♠♦♥❧② ✉s❡❞ t♦ ♠✐♥✐♠✐③❡ ❝♦♥str❛✐♥❡❞ ❢✉♥❝t✐♦♥s ◆♦t❡ ❛❧s♦ t❤❛t ✇❡ ❞♦♥✬t ♣❡♥❛❧✐③❡ t❤❡ ✐♥t❡r❝❡♣t✱ β ✵ ✱ ♦r t❤❡ ❧✐♥❡❛r ❝♦❡✣❝✐❡♥t✱ β ✶ ✽✼ ✴ ✶✻✵

❘✐❞❣❡❞ P❲▲ ❘❡❣r❡ss✐♦♥ ❈❤♦♦s✐♥❣ ❉ = ■ K ②✐❡❧❞s t❤❡ r✐❞❣❡ ♣❡♥❛❧t② � � �� ✷ � n � K + λ ❜ T ❜ , Q ( λ ) = Y i − β ✵ + β ✶ x i + b k ( x i − c k ) + i = ✶ k = ✶ � n � K � � �� ✷ + λ ① T i β + ③ T b ✷ = Y i − i ❜ k i = ✶ k = ✶ ( ❨ − ❳ β + ❩❜ ) T ( ❨ − ❳ β + ❩❜ ) + λ ❜ T ❜ = ✇❤❡r❡ ❳ ✐s ❛♥ n × ✷ ♠❛tr✐① [ ✶ , ① ] ✱ β = ( β ✵ , β ✶ ) T ✱ ❩ ✐s ❛♥ n × K s♣❧✐♥❡ ❜❛s✐s ❞❡s✐❣♥ ♠❛tr✐① ✇✐t❤ ( i , j ) ✲t❤ ❡❧❡♠❡♥t ❡q✉❛❧ t♦ ( x i − c j ) + ✱ ❛♥❞ ❜ = ( b ✶ , . . . , b K ) T ✳ ◆♦t❡✱ t❤❡ ( i , j ) ✲t❤ ❡❧❡♠❡♥t ♦❢ ❩ ij = ✵ ✐❢ x i ≤ c j ❛♥❞ ❡q✉❛❧ t♦ x i − c j ✐❢ x i > c j ❢♦r i = ✶ , . . . , n ❛♥❞ j = ✶ , . . . , K ✳ ✽✽ ✴ ✶✻✵

▼✐①❡❞ ▼♦❞❡❧ ❙tr✉❝t✉r❡ Q ( λ ) ❤❛s t❤❡ ❢♦r♠ ♦❢ ❛ ♠✐①❡❞ ♠♦❞❡❧ ✇✐t❤ ✜①❡❞ ❡✛❡❝ts β ❛♥❞ r❛♥❞♦♠ ❡✛❡❝ts ❜ ❲❡ ✇✐❧❧ ❞✐s❝✉ss ♠✐①❡❞ ♠♦❞❡❧s s❤♦rt❧②✱ ❜✉t ❢♦r ♥♦✇ ♥♦t❡ t❤❛t t❤❡ r❛♥❞♦♠ ❡✛❡❝ts ❛r❡ ♥♦t s✉❜❥❡❝t✲s♣❡❝✐✜❝ ✭✐✳❡✳✱ ♥♦t ❜ i ✮ ■♥st❡❛❞✱ t❤❡② ❛r❡ s❤❛r❡❞ ❜② ❛❧❧ n s✉❜❥❡❝ts ◆❡✈❡rt❤❡❧❡ss✱ ✇❡ ❝❛♥ ✉s❡ ❛ ♠✐①❡❞ ♠♦❞❡❧ ❢r❛♠❡✇♦r❦ ✭❢r❡q✉❡♥t✐st ♦r ❇❛②❡s✐❛♥✮ t♦ ❡st✐♠❛t❡ β ✱ λ ✱ ❜ ❛♥❞ σ ✷ ✽✾ ✴ ✶✻✵

▼✐①❡❞ ▼♦❞❡❧ ❙tr✉❝t✉r❡ ■♥ ♣❛rt✐❝✉❧❛r✱ ✇❡ ❝❛♥ ✇r✐t❡ ♦✉r ♠♦❞❡❧ ❛s n × ✶ = ❳ ❨ p × ✶ + ❩ n × K ❜ K × ✶ + ❡ n × ✶ , n × p β ✇❤❡r❡ β = ( β ✵ , β ✶ ) T ✐s ❛ p × ✶ ( p = ✷ ) ✈❡❝t♦r ❢♦r ✜①❡❞ ❡✛❡❝ts✱ ❜ ∼ ◆ ( ✵ , σ ✷ b ■ K ) ✐s ❛ ✈❡❝t♦r ♦❢ r❛♥❞♦♠ ❡✛❡❝ts✱ ❡ ∼ ◆ ( ✵ , σ ✷ e ■ n ) ✱ ❛♥❞ σ ✷ b = σ ✷ e /λ ✳ ◆ ( ✵ , σ ✷ b ■ K ) ✐s ❛ ♣r✐♦r ❞✐str✐❜✉t✐♦♥ ❢♦r ❜ t❤❛t ✐♠♣♦s❡s ♣r✐♦r s❤r✐♥❦❛❣❡ t♦✇❛r❞ ③❡r♦✳ ❋♦r ✜①❡❞ σ ✷ e ✱ ❛s λ → ✵✱ σ ✷ b → ∞ ✳ ❚❤✉s✱ t❤❡r❡ ✐s ❤✐❣❤ ✈❛r✐❛❜✐❧✐t② ✐♥ t❤❡ ♣r✐♦r✱ ❛♥❞ ❧✐tt❧❡ s❤r✐♥❦❛❣❡✴s♠♦♦t❤✐♥❣ ❆s λ → ∞ ✱ σ ✷ b → ✵✳ ❍❡r❡✱ t❤❡r❡ ✐s ❧♦✇ ✈❛r✐❛❜✐❧✐t②✱ ❛♥❞ ❜ ✐s ♠♦r❡ t✐❣❤t❧② ❝❡♥t❡r❡❞ ❛❜♦✉t ③❡r♦ ✭✐♥❝r❡❛s❡❞ s❤r✐♥❦❛❣❡✮ ✾✵ ✴ ✶✻✵

❋r❡q✉❡♥t✐st ■♥❢❡r❡♥❝❡ ❢♦r β ❲❡ ❝❛♥ s❤♦✇ ✹ t❤❛t ✐❢ σ ✷ e ❛♥❞ σ ✷ b ❛r❡ ❦♥♦✇♥✱ t❤❡♥ t❤❡ ❇▲❯❊ ♦❢ β ✐s � � − ✶ ❳ Σ − ✶ ❨ , � ❳ T Σ − ✶ ❳ p × ✶ = β b ❩❩ T ✐s t❤❡ ♠❛r❣✐♥❛❧ ✈❛r✐❛♥❝❡ ♦❢ ❨ ❛❢t❡r n × n = σ ✷ e ■ n + σ ✷ ✇❤❡r❡ Σ ✐♥t❡❣r❛t✐♥❣ ♦✉t ❜ ✳ ❚❤✐s ✐s ♦❜t❛✐♥❡❞ ❜② ♠❛①✐♠✐③✐♥❣ t❤❡ ♠❛r❣✐♥❛❧ ❧✐❦❡❧✐❤♦♦❞ ♦❢ ❨ ❛❢t❡r ✐♥t❡❣r❛t✐♥❣ ♦✉t t❤❡ r❛♥❞♦♠ ❡✛❡❝ts ❜ ✳ ❈❧❡❛r❧②✱ � β ✐s t❤❡ ✇❡✐❣❤t❡❞ ❧❡❛st sq✉❛r❡s ❡st✐♠❛t❡ ♦❢ � β ✇❤❡r❡ ❈♦✈ ( ❨ | ❳ ) = Σ ✳ ✹ ❙❡❡ ❲❛❦❡✜❡❧❞✱ ♣✳ ✺✻✺✕✺✻✻✳ ✾✶ ✴ ✶✻✵

❋r❡q✉❡♥t✐st ■♥❢❡r❡♥❝❡ ❢♦r ❜ ▲✐❦❡✇✐s❡✱ ✐❢ t❤❡ ✈❛r✐❛♥❝❡ ❝♦♠♣♦♥❡♥ts ❛r❡ ❦♥♦✇♥✱ t❤❡ ❜❡st ❧✐♥❡❛r ✉♥❜✐❛s❡❞ ♣r❡❞✐❝t♦r ✭❇▲❯P✮ ♦❢ ❜ ✐s � � − ✶ ❩ T ❘ − ✶ ( ❨ − ❳ � ● − ✶ + ❩ T ❘ − ✶ ❩ ˜ ❜ = β ) K × ✶ b ❩ T Σ − ✶ ( ❨ − ❳ � σ ✷ = β ) , ✇❤❡r❡ ● = σ ✷ b ■ k ✱ ❘ = σ ✷ e ■ n ✱ ❛♥❞ t❤❡ ❧❛st ❧✐♥❡ ❢♦❧❧♦✇s ❢r♦♠ t❤❡ ❢❛❝t t❤❛t✱ ❢♦r ✐♥✈❡rt✐❜❧❡ ♠❛tr✐❝❡s ❆ ❛♥❞ ❈ ✱ � � − ✶ ( ❆ + ❇❈❉ ) − ✶ ❇❈ ❈ − ✶ + ❉❆ − ✶ ❇ ❆ − ✶ ❇ = ✾✷ ✴ ✶✻✵

❯♥❦♥♦✇♥ ❱❛r✐❛♥❝❡ ❈♦♠♣♦♥❡♥ts ❙✐♥❝❡ σ ✷ e ❛♥❞ σ ✷ b ❛r❡ ✉♥❦♥♦✇♥✱ ✐♥ ♣r❛❝t✐❝❡ ✇❡ ♣❧✉❣ ✐♥ t❤❡ ❘❊▼▲ b ❛♥❞ � ❡st✐♠❛t❡s ♦❢ σ ✷ e ✱ σ ✷ β t♦ ❞❡r✐✈❡ ❧❛r❣❡✲s❛♠♣❧❡ ❡st✐♠❛t♦rs ❛♥❞ ♣r❡❞✐❝t♦rs � � − ✶ ❳ T � − ✶ ❳ − ✶ ❨ � ❳ � = ❛♥❞ β Σ Σ − ✶ ( ❨ − ❳ � b ❩ T � � σ ✷ ❜ = � Σ β ) ❋✐♥❛❧❧②✱ t❤❡ ❘❊▼▲ ❡st✐♠❛t❡ ♦❢ t❤❡ s♠♦♦t❤✐♥❣ ♣❛r❛♠❡t❡r ✐s ˆ σ ✷ σ ✷ λ = � e / � b ✳ ✾✸ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ■♥ t❤❡ ❇❛②❡s✐❛♥ s❡tt✐♥❣✱ ✇❡ ✉s❡ t❤❡ s❛♠❡ ♠✐①❡❞ ♠♦❞❡❧ ❢r❛♠❡✇♦r❦ ✇✐t❤ ❛ ◆ K ( ✵ , σ ✷ b ■ K ) ♣r✐♦r ❢♦r ❜ ❇✉t ♥♦✇ ✇❡ ♣❧❛❝❡ ❛❞❞✐t✐♦♥❛❧ ❝♦♥❥✉❣❛t❡ ♣r✐♦rs ♦♥ β ✱ σ ✷ e ❛♥❞ σ ✷ b ❖r✱ ❡q✉✐✈❛❧❡♥t❧②✱ ✇❡ ♣❧❛❝❡ ♣r✐♦rs ♦♥ τ e = ✶ /σ ✷ e ❛♥❞ τ b = ✶ /σ ✷ b ❆ ❝♦♠♠♦♥ ❝❤♦✐❝❡ ✐s✿ ✶ β ∼ ◆ ✷ ( β ✵ , ❚ − ✶ ✵ ) ✷ π ( σ ✷ e ) ∝ ✶ /σ ✷ e ✭✐♠♣r♦♣❡r ♣r✐♦r✮ ✸ τ b ∼ ●❛ ( c , d ) ❚❤❡ ❝❤♦✐❝❡ ♦❢ c ❛♥❞ d ❛✛❡❝ts s♠♦♦t❤♥❡ss✳ ❚♦ s❡❧❡❝t ✈❛❧✉❡s✱ ❝❛♥ ✉s❡ ❇❛②❡s✐❛♥ ✐♥❢♦r♠❛t✐♦♥ ❝r✐t❡r✐❛✱ ♦r ♣❧❛❝❡ ❛ ❞✐s❝r❡t❡ ❜✐✈❛r✐❛t❡ ♣r✐♦r ♦♥ ( c , d ) ✾✹ ✴ ✶✻✵

❇❛②❡s✐❛♥ ■♥❢❡r❡♥❝❡ ∗ ❚❤✐s ❧❡❛❞s t♦ t❤❡ ❢♦❧❧♦✇✐♥❣ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧s✿ ✶ β | ② , r❡st ∼ ◆ ✷ ( ♠ , ❱ ) ✱ ✇❤❡r❡ � � − ✶ ✱ ✇❤❡r❡ τ e = ✶ /σ ✷ ❚ ✵ β ✵ + τ e ❳ T ❳ ✷ × ✷ = ❱ e � � ❚ ✵ β ✵ + τ e ❳ T ( ② − ❩❜ ) ✷ × ✶ = ❱ ♠ ✷ ❜ | ② , r❡st ∼ ◆ K ( ♠ , ❱ ) ✱ ✇❤❡r❡ � � − ✶ τ b ■ K + τ e ❩ T ❩ K × K = ❱ K × ✶ = τ e ❱ ❩ T ( ② − ❳ β ) ♠ � n � ✷ , ( ② − ❳ β − ❩❜ ) T ( ② − ❳ β − ❩❜ ) ✸ τ e | ② , r❡st ∼ ●❛ ✷ � � c + K ✷ , d + ❜ T ❜ / ✷ ✹ τ b | ② , r❡st ∼ ●❛ ❙❡❡ ♣❡♣t✐❞❡✳r ❢♦r ❞❡t❛✐❧s ✾✺ ✴ ✶✻✵

❘ ❈♦❞❡ ❢♦r P❡♥❛❧✐③❡❞ P❲▲ ▼♦❞❡❧ Gibbs Sampler for Penalized PWL Spline Model # Import Data: x=age, y=logc # PWL splines k<-10 # number of knots grid<- seq (5,14,length=k) # k+1 initial knot locations (includes boundaries) Z<- matrix (0,n,k) for (j in 1:k) Z[,j]<- pmax (age-grid[j],0) # Priors c<-d<-1e-5 # Hyperpriors for taub -- increase to 1 to reduce smoothness # Prior var of taub small and var sigma_b is large = less shrinkage T0<- diag (.0001,p) # Prior Precision # Inits -- see R code # Fine grid of x values with which to plot yhat x2<- seq (0,16,by=.01) # Grid of age values X2<- cbind (1,x2) # Fixed effect covs Z2<- matrix (0, length (x2),k) for (j in 1:k) Z2[,j]<- pmax (x2-grid[j],0) ######################### # Gibbs of Ridged Model # ######################### for (i in 1:nsim){ # Update beta v<- solve (T0+taue* crossprod (X)) m<-v%*%(T0%*%beta0+taue* t (X)%*%(y-Z%*%b)) Beta[i,]<-beta<- c ( rmvnorm (1,m,v)) # Update b v<- solve (taub* diag (k)+taue* crossprod (Z)) m<-taue*v%*% t (Z)%*%(y-X%*%beta) B[i,]<-b<- c ( rmvnorm (1,m,v)) # Update taue (diffuse) g<- crossprod (y-X%*%beta-Z%*%b) taue<- rgamma (1,n/2,g/2) Sigmae[i]<-1/taue # Update taub (proper) taub<- rgamma (1,c+k/2,d+ crossprod (b)/2) Sigmab[i]<-1/taub # Yhat Yhat[i,]<-X2%*%beta+Z2%*%b # Predicted values } ✾✻ ✴ ✶✻✵

❘❊▼▲ ❊st✐♠❛t❡ ♦❢ ▼❡❛♥ ❘❡❣r❡ss✐♦♥ ❋✉♥❝t✐♦♥ ❋✐❣✉r❡ ✹✿ ❘❊▼▲ ❡st✐♠❛t❡s ♦❢ t❤❡ ♦r✐❣✐♥❛❧ ❛♥❞ r✐❞❣❡✲s♠♦♦t❤❡❞ P❲▲ r❡❣r❡ss✐♦♥ ❢✉♥❝t✐♦♥s✳ 0.8 Log Concentration 0.7 0.6 0.5 0 5 10 15 Age (Years) Log Concentration Orginal PWL Ridged PWL Lower CI Ridged PWL Upper CI Ridged PWL ✾✼ ✴ ✶✻✵

P♦st❡r✐♦r ▼❡❛♥ ❊st✐♠❛t❡ ♦❢ t❤❡ ❘❡❣r❡ss✐♦♥ ❋✉♥❝t✐♦♥ ❋✐❣✉r❡ ✺✿ P♦st❡r✐♦r ♠❡❛♥ ❡st✐♠❛t❡s ♦❢ t❤❡ ♦r✐❣✐♥❛❧ ❛♥❞ r✐❞❣❡✲s♠♦♦t❤❡❞ P❲▲ r❡❣r❡ss✐♦♥ ❢✉♥❝t✐♦♥s✳ Original PWL ● 0.80 Ridged 95% CrI ● 0.75 ● ● ● ● ● ● ● ● 0.70 ● ● ● ●● ● ● Log Concentration ● ● ● ● ● ● ● ● ● ● ● ● 0.65 ● ● ● ● ● 0.60 ● ● ● ● ● 0.55 ● ● 0.50 ● 0 5 10 15 Age (Y ears) ✾✽ ✴ ✶✻✵

P♦st❡r✐♦r ▼❡❛♥ ❊st✐♠❛t❡s ✇✐t❤ ❤②♣❡r♣❛r❛♠❡t❡rs c = d = ✶ Original PWL 0.80 Ridged 95% CrI 0.75 0.70 Log Concentration 0.65 0.60 0.55 0.50 0 5 10 15 Age (Y ears) ✾✾ ✴ ✶✻✵

▲✐♥❡❛r ▼✐①❡❞ ▼♦❞❡❧s ❈♦♥s✐❞❡r ❛ ❧✐♥❡❛r ♠✐①❡❞ ♠♦❞❡❧ ♦❢ t❤❡ ❢♦r♠ n i × ✶ = ❳ i ❨ i n i × p β p × ✶ + ❩ i n i × q ❜ i q × ✶ + ❡ i n i × ✶ , i = ✶ , . . . , n ✇❤❡r❡✱ ✉♥❞❡r ❝❧❛ss✐❝❛❧ ❛ss✉♠♣t✐♦♥s✱ • ❜ i ∼ ◆ q ( ✵ , ● ) ✐s ❛ ✈❡❝t♦r ♦❢ s✉❜❥❡❝t✲s♣❡❝✐✜❝ r❛♥❞♦♠ ❡✛❡❝ts ✭❝♦✉❧❞ ❛❧❧♦✇ ❢♦r ♥♦♥✲♥♦r♠❛❧✐t②✮ • ❩ i ✐s ❛ r❛♥❞♦♠ ❡✛❡❝t ❞❡s✐❣♥ ♠❛tr✐① • ❡ i ∼ ◆ n i ( ✵ , ❘ i ) ∼ ◆ n i ( ✵ , σ ✷ e ■ n i ) • ❜ i ❛♥❞ ❡ i ❛r❡ ✐♥❞❡♣❡♥❞❡♥t ❘❱s ◆♦t❡ t❤❛t ✇❡ ❛❧❧♦✇ ❢♦r ✉♥❜❛❧❛♥❝❡❞ ❞❛t❛ ✇✐t❤ t❤❡ ♥✉♠❜❡r ♦❢ r❡♣❡❛t❡❞ ♠❡❛s✉r❡♠❡♥ts✱ n i ✱ ✈❛r②✐♥❣ ❛❝r♦ss s✉❜❥❡❝ts ✶✵✵ ✴ ✶✻✵