s trts r - - PowerPoint PPT Presentation
s trts r r r s Prt tr rt
✶ ✴ ✶✻✵
i β + ei i = ✶, . . . , n,
✐✐❞
e )
e ✐s ❛ ♣r❡❝✐s✐♦♥ t❡r♠
✷ ✴ ✶✻✵
e
e
✸ ✴ ✶✻✵
e ✭♦r✱
✵ ), ✇❤❡r❡ ❚ ✵ ✐s t❤❡ ♣r✐♦r ♣r❡❝✐s✐♦♥
e ✐s ✐♥✈❡rs❡✲●❛♠♠❛ ✭■●✮ ✇✐t❤
✹ ✴ ✶✻✵
✵
✺ ✴ ✶✻✵
✻ ✴ ✶✻✵
✼ ✴ ✶✻✵
library(mvtnorm) # Priors beta0<-rep(0,p) # Prior mean of beta, where p=# of parameters T0<-diag(0.01,p) # Prior Precision of beta a<-b<-0.001 # Gamma hyperparms for taue # Inits taue<-1 # Error precision # MCMC Info nsim<-1000 # Number of MCMC Iterations thin<-1 # Thinning interval burn<-nsim/2 # Burnin lastit<-(nsim-burn)/thin # Last stored value # Store Beta<-matrix(0,lastit,p) # Matrices to store results Sigma2<-rep(0,lastit) Resid<-matrix(0,lastit,n) # Store resids Dy<-matrix(0,lastit,512) # Store density values for residual density plot Qy<-matrix(0,lastit,100) # Store quantiles for QQ plot ######### # Gibbs # ######### tmp<-proc.time() # Store current time for (i in 1:nsim){ # Update beta v<-solve(T0+taue*crossprod(X)) m<-v%*%(T0%*%beta0+taue*crossprod(X,y)) beta<-c(rmvnorm(1,m,v)) # Update tau taue<-rgamma(1,a+n/2,b+crossprod(y-X%*%beta)/2) # Store Samples if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Sigma2[j]<-1/taue Resid[j,]<-resid<-y-X%*%beta # Raw Resid Dy[j,]<-density(resid/sd(resid),from=-5,to=5)$y # Density of Standardized Resids Qy[j,]<-quantile(resid/sd(resid),probs=seq(.001,.999,length=100)) # Quantiles for QQ Plot } if (i%%100==0) print(i) } run.time<-proc.time()-tmp # MCMC run time # Took 1 second to run 1000 iterations with n=1000 subjects
✽ ✴ ✶✻✵
iid
e)
e
✾ ✴ ✶✻✵
−4 −2 2 4 0.0 0.1 0.2 0.3 0.4
Density Plot of Standardized Residuals
Quantile Density MLE Bayes
−2 −1 1 2 3 −3 −1 1 2 3
QQ Plot of Standardized Residuals
Normal Quantile Sample Quantiles
Bayes ✶✵ ✴ ✶✻✵
✶❖✬❍❛❣❛♥ ❛♥❞ ▲❡♦♥❛r❞ ✭✶✾✼✻✮❀ ❆③③❛❧✐♥✐✱ ✶✾✽✺
✶✶ ✴ ✶✻✵
−15 −10 −5 0.00 0.05 0.10 0.15
True Errors (e)
y Density α = −5 ω = 4 −4 −2 2 0.0 0.1 0.2 0.3 0.4
Density of Standardized Residuals
Standardized Residual Density MLE Bayes
−2 −1 1 2 3 −3 −2 −1 1 2 3
QQ Plot of Standardized Residuals
Normal Quantile Observed Quantile
Bayes
✶✷ ✴ ✶✻✵
i β,
i β, ✶) ❛♥❞
i β)
✶✸ ✴ ✶✻✵
−2 2 4 0.0 0.1 0.2 0.3 0.4 zi f(zi) Pr(Zi > 0) = Pr(Yi = 1)
xi
Tβ ✶✹ ✴ ✶✻✵
✶✺ ✴ ✶✻✵
✶ ❉r❛✇ ③ ❢r♦♠ π(③|②, β) ✷ ❉r❛✇ β ❢r♦♠ π(β|③)
✶✻ ✴ ✶✻✵
✵
✶ ❋♦r i = ✶, . . . , n✱ s❛♠♣❧❡ zi ❢r♦♠ ❛ ◆(①T i β, ✶) ❞✐str✐❜✉t✐♦♥
✷ ❙❛♠♣❧❡ β ❢r♦♠ ◆p(♠, ❱ )✱ ✇❤❡r❡ ♠ = ❱
✶✼ ✴ ✶✻✵
# 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
✶✽ ✴ ✶✻✵
† ❇❛s❡❞ ♦♥ ❆❧❜❡rt ❛♥❞ ❈❤✐❜ ❞❛t❛ ❛✉❣♠❡♥t❛t✐♦♥ s❛♠♣❧❡r✳
✶✾ ✴ ✶✻✵
✷✵ ✴ ✶✻✵
zi f(zi) Pr(Zi > 0) = Pr(Yi = 1)
xi
Tβ t3(xi
Tβ)
N(xi
Tβ,1)
0.0 0.1 0.2 0.3 0.4 ✷✶ ✴ ✶✻✵
−4 −2 2 4 −5 5 t quantile logistic quantile t3 t8 t16 Logistic*0.634
✷✷ ✴ ✶✻✵
✶ ●❡♥❡r❛t✐♥❣ λi ❢r♦♠ ●❛(ν/✷, ν/✷) ✭♣❛r❛♠❡t❡r✐③❡❞ ❛s r❛t❡✮ ✷ ●❡♥❡r❛t✐♥❣ Zi|λi ❢r♦♠ ◆
i β, λ−✶ i
i β) ❛♥❞ t❤❡ ♦r✐❣✐♥❛❧ ❜✐♥❛r②
✷✸ ✴ ✶✻✵
✶ ❋♦r i = ✶, . . . , n✱ s❛♠♣❧❡ zi ❢r♦♠ ✐ts tr✉♥❝❛t❡❞ ◆(①T i β, λ−✶ i ) ❞✐str✐❜✉t✐♦♥
✷ ❙❛♠♣❧❡ β ❢r♦♠ ◆p(♠, ❱ )✱ ✇❤❡r❡
✸ ❋♦r i = ✶, . . . , n✱ s❛♠♣❧❡ λi ❢r♦♠ ●❛
i β)✷)/✷
✺ ❋♦r t❤❡ ❧♦❣✐st✐❝ ❛♣♣r♦①✐♠❛t✐♦♥✱ s❡t ν = ✽ ❛♥❞ ❞✐✈✐❞❡ t❤❡ ♣♦st❡r✐♦r ♠❡❛♥s ❛♥❞ ❙❉s
✷✹ ✴ ✶✻✵
# 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)
✷✺ ✴ ✶✻✵
† ❇❛s❡❞ ♦♥ ❆❧❜❡rt ❛♥❞ ❈❤✐❜ t✲❧✐♥❦ ♠♦❞❡❧✳
✷✻ ✴ ✶✻✵
d
∞
✷✼ ✴ ✶✻✵
0.0 0.5 1.0 1.5 1 2 3 4 x f(x) PG(1, 0) Ga(2, 10) ✷✽ ✴ ✶✻✵
✵
✷✾ ✴ ✶✻✵
n
n
n
i β✳
✸✵ ✴ ✶✻✵
✵
i /✷p(ωi|✶, ✵) ❞ωi,
✷✳
i β✱ ❛♥❞ t❤❡♥ ❛✈❡r❛❣❡ ❛❝r♦ss t❤❡
✸✶ ✴ ✶✻✵
ωi
✸✷ ✴ ✶✻✵
✵ ) ♣r✐♦r ❢♦r β✱ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧
✶ ❋♦r i = ✶, . . . , n✱ ✉♣❞❛t❡ ωi ❢r♦♠ ❛ P●(✶, ηi) ❞❡♥s✐t②✱ ✇❤❡r❡
i β ✷ ❋♦r i = ✶, . . . , n✱ ❞❡✜♥❡ zi = yi−✶/✷ ωi ✸ ❈♦♥❞✐t✐♦♥❛❧ ♦♥ ③ ❛♥❞ ❲ ✱ ✉♣❞❛t❡ β ❢r♦♠ ◆p(♠, ❱ )✱ ✇❤❡r❡
✸✸ ✴ ✶✻✵
✸✹ ✴ ✶✻✵
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
✸✺ ✴ ✶✻✵
† ❇❛s❡❞ ♦♥ P● s❛♠♣❧❡r✳
✸✻ ✴ ✶✻✵
i β,
✸✼ ✴ ✶✻✵
zi f(zi) 0.0 0.1 0.2 0.3 0.4
Pr(Zi < γ1) = Pr(Yi = 1) Pr(γ1 < Zi<γ2) = Pr(Yi = 2) Pr(Zi > γ2) = Pr(Yi = 3)
xi
Tβ
γ1 γ2
✸✽ ✴ ✶✻✵
✸✾ ✴ ✶✻✵
✶ ❋♦r ✐❂✶✱✳ ✳ ✳ ✱♥✱ ✐❢ yi = k✱ ❞r❛✇ zi ❢r♦♠ ◆(①T i β, λ−✶ i ) tr✉♥❝❛t❡❞ ❜❡t✇❡❡♥ γk−✶ ❛♥❞
✷ ❯♣❞❛t❡ β ❢r♦♠ ◆p(♠, ❱ )✱
✸ ❢♦r k = ✶, . . . , K − ✶✱ ✉♣❞❛t❡ γk ❢r♦♠ ❯♥✐❢(♠❛① zi : yi = k, ♠✐♥ zi : yi = k + ✶) ✹ ❋♦r t❤❡ ❝✉♠✉❧❛t✐✈❡ ❧♦❣✐t ♠♦❞❡❧✱ ✉♣❞❛t❡ λi ❢r♦♠ ✐ts ●❛♠♠❛ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧
✹✵ ✴ ✶✻✵
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
✹✶ ✴ ✶✻✵
† ❇❛s❡❞ ♦♥ ❆❧❜❡rt ❛♥❞ ❈❤✐❜ s❛♠♣❧❡r✳
✹✷ ✴ ✶✻✵
100 200 300 400 500 −1.4 −1.0 −0.6 Iteration γ1 100 200 300 400 500 0.6 0.8 1.0 1.2 Iteration γ2 100 200 300 400 500 0.6 0.8 Iteration β
✹✸ ✴ ✶✻✵
k=✶ πik = ✶
i βk
j=✶ ❡①T
i βj , k = ✶, . . . , K,
i βk
j=✶ ❡①T
i βj , k = ✶, . . . , K − ✶
j=✶ ❡①T
i βj
k=✶ πik = ✶✳ ✹✹ ✴ ✶✻✵
i βk),
i βk.
✹✺ ✴ ✶✻✵
✹✻ ✴ ✶✻✵
✹✼ ✴ ✶✻✵
i βk
j=✶ ❡①T
i βj , k = ✶, . . . , K,
✹✽ ✴ ✶✻✵
n
ik (✶ − πik)✶−Uik
❡①T
i βk
K
j=✶ ❡①T i βj ✹✾ ✴ ✶✻✵
i βk−cik
i βk−cik
j=k ❡①T
i βj ❛♥❞ ηik = ①T
i βk − cik✳
j=k ❡①T
i βj ✐♥❝❧✉❞❡s t❤❡ r❡❢❡r❡♥❝❡ ❝❛t❡❣♦r② K✳
i βK = ✶ ❛♥❞ ❤❡♥❝❡
i βj = ❧♦❣
∈{k, K}
i βj
✺✵ ✴ ✶✻✵
n
n
✺✶ ✴ ✶✻✵
✵ ) ❢♦r β✶, . . . , βK−✶✱ t❤❡ ●✐❜❜s s❛♠♣❧❡r ♣r♦❝❡❡❞s
✶ ❖✉ts✐❞❡ ♦❢ t❤❡ ●✐❜❜s ❧♦♦♣✱ ❞❡✜♥❡ uik = ✶(yi=k) ❢♦r i = ✶, . . . , n ❛♥❞
✷ ❋♦r i = ✶, . . . , n ❛♥❞ k = ✶, . . . , K − ✶✱ ✉♣❞❛t❡ ωik ❢r♦♠ ❛
i βk − cik✱ ❛♥❞ cik ✇❛s ❞❡✜♥❡❞
✸ ❋♦r i = ✶, . . . , n ❛♥❞ k = ✶, . . . , K − ✶✱ ❞❡✜♥❡ zik = uik−✶/✷ ωik
✹ ❋♦r k = ✶, . . . , K − ✶✱ ✉♣❞❛t❡ βk ❢r♦♠ ◆p(♠k, ❱ k)✱ ✇❤❡r❡
✺✷ ✴ ✶✻✵
i βk
j=✶ ❡①T
i βj , k = ✶, . . . , ✸,
✺✸ ✴ ✶✻✵
# 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)
✺✹ ✴ ✶✻✵
# 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❡❞ ♦♥ P♦❧s♦♥ ❡t ❛❧✳ Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r✳
✺✻ ✴ ✶✻✵
d
i , r > ✵,
i β
i β
✺✼ ✴ ✶✻✵
✺✽ ✴ ✶✻✵
✵
i /✷p(ωi|r + yi, ✵) ❞ωi,
✺✾ ✴ ✶✻✵
✻✵ ✴ ✶✻✵
✶ ❋♦r i = ✶, . . . , n✱ ❞r❛✇ ωi ❢r♦♠ ✐ts P●(yi + r, ηi) ❞✐str✐❜✉t✐♦♥✱
i β ✷ ❋♦r i = ✶, . . . , n✱ ❞❡✜♥❡ zi = yi − r
✸ ❆ss✉♠✐♥❣ ❛ ◆p
✵
✻✶ ✴ ✶✻✵
❩❤♦✉ ❛♥❞ ❈❛r✐♥ ✭✷✵✶✺✮ ❛♥❞ ❉❛❞❛♥❡❤ ❡t ❛❧✳ ✭✷✵✶✽✮ ❞❡s❝r✐❜❡ ❛ t✇♦✲st❡♣ ❝♦♥❥✉❣❛t❡ ●✐❜❜s ✉♣❞❛t❡ ❢♦r r ❚❤❡ ❛♣♣r♦❛❝❤ ✐♥tr♦❞✉❝❡s ❛ s❛♠♣❧❡ ♦❢ ❧❛t❡♥t ❝♦✉♥ts✱ li✱ ✉♥❞❡r❧②✐♥❣ ❡❛❝❤ ♦❜s❡r✈❡❞ ❝♦✉♥t yi ❈♦♥❞✐t✐♦♥❛❧ ♦♥ yi ❛♥❞ r✱ li ❤❛s ❛ ❞✐str✐❜✉t✐♦♥ ❞❡✜♥❡❞ ❜② ❛ ❈❤✐♥❡s❡ r❡st❛✉r❛♥t t❛❜❧❡ ✭❈❘❚✮ ❞✐str✐❜✉t✐♦♥✿ li =
yi
uj uj ∼ ❇❡r♥
r + j − ✶
❚❤❡ ♥❛♠❡ ✏❈❤✐♥❡s❡ r❡st❛✉r❛♥t t❛❜❧❡✑ ❞❡r✐✈❡s ❢r♦♠ t❤❡ ❢❛❝t t❤❛t uj = ✶ ✐❢ ❛ ♥❡✇ ❝✉st♦♠❡r s✐ts ❛♥❞ ❛♥ ✉♥♦❝❝✉♣✐❡❞ t❛❜❧❡ ✐♥ ❛ ❈❤✐♥❡s❡ r❡st❛✉r❛♥t ✭❛❝❝♦r❞✐♥❣ t♦ ❛ s♦✲❝❛❧❧❡❞ ✏❈❤✐♥❡s❡ r❡st❛✉r❛♥t ♣r♦❝❡ss✑✮✱ ❛♥❞ li ✐s t❤❡ t♦t❛❧ ♥✉♠❜❡r ♦❢ ♦❝❝✉♣✐❡❞ t❛❜❧❡s ✐♥ t❤❡ r❡st❛✉r❛♥t ❛❢t❡r yi ❝✉st♦♠❡rs ❙♦ ✐♥ ❙t❡♣ ✶✱ ✇❡ ❞r❛✇ li (i = ✶, . . . , n) ❛❝❝♦r❞✐♥❣ t♦ t❤✐s ❈❘❚ ❞✐str✐❜✉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✱ li
ind
∼ P♦✐[−r ❧♥(✶ − ψi)], ✇❤❡r❡ ψi = ❡①♣
i β
i β
, i = ✶, . . . , n. ❙❡❡ ❉❛❞❛♥❡❤ ❡t ❛❧✳ ✭✷✵✶✽✮ ❛♥❞ ❩❤♦✉ ❛♥❞ ❈❛r✐♥ ✭✷✵✶✺✮ ❢♦r ❞❡t❛✐❧s ❚❤✉s✱ ✐❢ ✇❡ ❛ss✉♠❡ ❛ ●❛(a, b) ♣r✐♦r ❢♦r r✱ t❤❡♥ t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r r ✐♥ ❙t❡♣ ✷ ✐s r|❧, ψ ∼
n
li, b −
n
❧♥(✶ − ψi)
✇❤❡r❡ ❧ = (l✶, . . . , ln)T ❛♥❞ ψ = (ψ✶, . . . , ψn)T✳ ❚❤❡ ●✐❜❜s ✉♣❞❛t❡ ✜rst ❞r❛✇s li (i = ✶, . . . , n) ✐♥❞❡♣❡♥❞❡♥t❧② ❢r♦♠ ❛ ❈❘❚ ❞✐str✐❜✉t✐♦♥✱ ❛♥❞ t❤❡♥ r ❢r♦♠ ✐ts ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ●❛♠♠❛ ❞✐str✐❜✉t✐♦♥ ❣✐✈❡♥ ❧ ❛♥❞ ψ
✻✸ ✴ ✶✻✵
# Priors (diffuse prior for r) beta0<-rep(0,p) T0<-diag(.01,p) s<-0.01 # Proposal variance
# 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 }
✻✹ ✴ ✶✻✵
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 }
✻✺ ✴ ✶✻✵
d
i , r > ✵,
† Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ▼❍ ✉♣❞❛t❡ ❢♦r r✳ ❆❝❝❡♣t❛♥❝❡ r❛t❡ ❂ ✸✼✪✳ ‡ Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ●✐❜❜s ✉♣❞❛t❡ ❢♦r r✳ ✻✻ ✴ ✶✻✵
200 400 600 800 1000 1.6 1.8 2.0 2.2 Iteration β1 200 400 600 800 1000 0.2 0.4 0.6 0.8 Iteration β2 200 400 600 800 1000 0.8 1.0 1.2 1.4 Iteration r
✻✼ ✴ ✶✻✵
✻✽ ✴ ✶✻✵
i , yi = ✶, ✷, . . .
i β
i β
✻✾ ✴ ✶✻✵
✼✵ ✴ ✶✻✵
i α = η✶i
d
i
i β
i β
✼✶ ✴ ✶✻✵
✶ ●✐✈❡♥ ❝✉rr❡♥t ♣❛r❛♠❡t❡r ✈❛❧✉❡s✱ ✉♣❞❛t❡ α ✉s✐♥❣ t❤❡ ●✐❜❜s
✷ ❈♦♥❞✐t✐♦♥❛❧ ♦♥ Wi = ✶✱ ✉♣❞❛t❡ β ✉s✐♥❣ t❤❡ ◆❇ ●✐❜❜s
✸ ❯♣❞❛t❡ r ✉s✐♥❣ ❛ r❛♥❞♦♠✲✇❛❧❦ ▼❡tr♦♣♦❧✐s✲❍❛st✐♥❣s st❡♣ ♦r
✹ ❯♣❞❛t❡ t❤❡ ❧❛t❡♥t ❛t✲r✐s❦ ✐♥❞✐❝❛t♦rs✱ W✶, . . . , Wn✱ ❢r♦♠ t❤❡✐r
✼✷ ✴ ✶✻✵
✼✸ ✴ ✶✻✵
i
i ),
✼✹ ✴ ✶✻✵
2 4 6 8 11 14 17 20 23 26 29 32 35 38 44 60 Count Percent 10 20 30 40 50 60 AIC for ZINB is 104 points lower than for NB 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 }
✼✻ ✴ ✶✻✵
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 }
✼✼ ✴ ✶✻✵
d
i
❚❛❜❧❡ ✽✿ ❘❡s✉❧ts ❢♦r ❩■◆❇ ▼♦❞❡❧✳
† Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ▼❍ ✉♣❞❛t❡ ❢♦r r✳ ❆❝❝❡♣t❛♥❝❡ r❛t❡ ❂ ✹✵✪✳ ‡ Pó❧②❛✲●❛♠♠❛ s❛♠♣❧❡r ✇✐t❤ ●✐❜❜s ✉♣❞❛t❡ ❢♦r r✳
✼✽ ✴ ✶✻✵
200 400 600 800 1000 −0.9 −0.7 −0.5 −0.3 1:lastit α1 200 400 600 800 1000 0.2 0.6 1.0 1:lastit α2 200 400 600 800 1000 1.6 1.8 2.0 2.2 1:lastit β1 200 400 600 800 1000 0.2 0.4 0.6 0.8 1:lastit β2 200 400 600 800 1000 0.8 1.0 1.2 1.4 1:lastit r
✼✾ ✴ ✶✻✵
200 400 600 800 1000 −1.0 −0.6 −0.2 1:lastit α1 200 400 600 800 1000 0.0 0.4 0.8 1:lastit α2 200 400 600 800 1000 1.6 1.8 2.0 2.2 1:lastit β1 200 400 600 800 1000 0.2 0.4 0.6 0.8 1:lastit β2 200 400 600 800 1000 0.8 1.0 1.2 1.4 1:lastit r
✽✵ ✴ ✶✻✵
✽✶ ✴ ✶✻✵
❋✐❣✉r❡ ✷✿ ❙❝❛tt❡r♣❧♦t ♦❢ ❧♦❣ ❈✲♣❡♣t✐❞❡ ❝♦♥❝❡♥tr❛t✐♦♥ ❜② ❛❣❡✳
ears) Log Concentration 2 4 6 8 10 12 14 16 0.5 0.6 0.7 0.8
✷❝✳❢✳✱ ❋✐t③♠❛✉r✐❝❡ ❡t ❛❧✳✱ ❈❤❛♣t❡r ✶✾
✽✷ ✴ ✶✻✵
K
✽✸ ✴ ✶✻✵
ears) Log Concentration 2 4 6 8 10 12 14 16 0.5 0.6 0.7 0.8
✽✹ ✴ ✶✻✵
✽✺ ✴ ✶✻✵
n
K
✶ ♠❛① |bk| < t ✷ K k=✶ |bk| < t ✭▲❛ss♦ ❝♦♥str❛✐♥t✮ ✸ K k=✶ b✷ k < t = ❜T❜ < t ✭❘✐❞❣❡ ❝♦♥str❛✐♥t✮ ✹ ❜T❉❜ < t ❢♦r s♦♠❡ K × K ♣♦s✐t✐✈❡ s❡♠✐✲❞❡✜♥✐t❡ ♣❡♥❛❧t②
✸❘✉♣♣❡rt ❡t ❛❧✳✱ ♣❛❣❡ ✻✺✳
✽✻ ✴ ✶✻✵
n
K
✽✼ ✴ ✶✻✵
n
K
n
i β + ③T i ❜
K
k
✽✽ ✴ ✶✻✵
✽✾ ✴ ✶✻✵
n×✶ = ❳ n×p β p×✶ + ❩ n×K ❜ K×✶ + ❡ n×✶,
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
e✱ ❛s λ → ✵✱ σ✷ b → ∞✳ ❚❤✉s✱ t❤❡r❡ ✐s ❤✐❣❤ ✈❛r✐❛❜✐❧✐t②
b → ✵✳ ❍❡r❡✱ t❤❡r❡ ✐s ❧♦✇ ✈❛r✐❛❜✐❧✐t②✱ ❛♥❞ ❜ ✐s ♠♦r❡
✾✵ ✴ ✶✻✵
e ❛♥❞ σ✷ b ❛r❡ ❦♥♦✇♥✱ t❤❡♥ t❤❡ ❇▲❯❊ ♦❢ β
p×✶ =
n×n = σ✷ e■ n + σ✷ b❩❩ T ✐s t❤❡ ♠❛r❣✐♥❛❧ ✈❛r✐❛♥❝❡ ♦❢ ❨ ❛❢t❡r
✹❙❡❡ ❲❛❦❡✜❡❧❞✱ ♣✳ ✺✻✺✕✺✻✻✳
✾✶ ✴ ✶✻✵
K×✶
b❩ TΣ−✶(❨ − ❳
b■ k✱ ❘ = σ✷ e■ n✱ ❛♥❞ t❤❡ ❧❛st ❧✐♥❡ ❢♦❧❧♦✇s ❢r♦♠ t❤❡
✾✷ ✴ ✶✻✵
e ❛♥❞ σ✷ b ❛r❡ ✉♥❦♥♦✇♥✱ ✐♥ ♣r❛❝t✐❝❡ ✇❡ ♣❧✉❣ ✐♥ t❤❡ ❘❊▼▲
e✱ σ✷ b ❛♥❞
−✶❳
−✶❨
b❩ T
−✶(❨ − ❳
e/
b✳
✾✸ ✴ ✶✻✵
b■ K) ♣r✐♦r ❢♦r ❜
e ❛♥❞ σ✷ b
e ❛♥❞ τb = ✶/σ✷ b
✶ β ∼ ◆✷(β✵, ❚ −✶ ✵ ) ✷ π(σ✷ e) ∝ ✶/σ✷ e ✭✐♠♣r♦♣❡r ♣r✐♦r✮ ✸ τb ∼ ●❛(c, d)
✾✹ ✴ ✶✻✵
✶ β|②, r❡st ∼ ◆✷(♠, ❱ )✱ ✇❤❡r❡
❱
✷×✷ =
−✶✱ ✇❤❡r❡ τe = ✶/σ✷
e
♠
✷×✶ = ❱
❱
K×K =
−✶ ♠
K×✶ = τe❱ ❩ T (② − ❳β) ✸ τe|②, r❡st ∼ ●❛
✾✺ ✴ ✶✻✵
# 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 }
✾✻ ✴ ✶✻✵
5 10 15
Age (Years)
0.5 0.6 0.7 0.8
Log Concentration
Upper CI Ridged PWL Lower CI Ridged PWL Ridged PWL Orginal PWL Log Concentration
✾✼ ✴ ✶✻✵
10 15 0.50 0.55 0.60 0.65 0.70 0.75 0.80 Age (Y ears) Log Concentration Original PWL Ridged 95% CrI ✾✽ ✴ ✶✻✵
5 10 15 0.50 0.55 0.60 0.65 0.70 0.75 0.80 Age (Y ears) Log Concentration Original PWL Ridged 95% CrI
✾✾ ✴ ✶✻✵
ni×✶ = ❳ i ni×p β p×✶ + ❩ i ni×q ❜i q×✶ + ❡i ni×✶,
e■ ni)
✶✵✵ ✴ ✶✻✵
e) = ◆ni(❳ iβ + ❩ i❜i, ❘i).
e) = ◆ni(❳ iβ, Σi),
i = σ✷ e■ ni + ❩ i●❩ T i ✳
✶✵✶ ✴ ✶✻✵
p×✶ =
i
−✶ i ❳ i
n
i
e ❛♥❞
p×p
i
−✶ i ❳
✶✵✷ ✴ ✶✻✵
i
✶✵✸ ✴ ✶✻✵
e) = ❊(˜
e)
e ❛r❡ ❦♥♦✇♥
✶✵✹ ✴ ✶✻✵
−✶ i
−✶ i
−✶ i ✱
✶✵✺ ✴ ✶✻✵
i Pi❩ i●,
i (■ ni − ❍i)
i Σ−✶ i ❳ i
i Σ−✶ i
n×✶) = σ✷ e(■ n − ❍)✳
✶✵✻ ✴ ✶✻✵
✶✵✼ ✴ ✶✻✵
i Pi❩ i●,
✶✵✽ ✴ ✶✻✵
e, ●) = ◆q(˜
i Pi❩ i●),
i (■ ni − ❍i)✳
✶✵✾ ✴ ✶✻✵
✶
e, ●)
i Σ−✶ i ❩ i●
i ❘−✶ i ❩ i
i ❩ i
e
❜i ❣✐✈❡♥ ② i✱ β✱ Σi ❛♥❞ ●
♦❢ ❛ ●✐❜❜s s❛♠♣❧❡r
✷ ❱❛r(˜
i Σ−✶ i (❨ i − ❳ i
i Pi❩ i●
β
♦❢ t❤❡ ♠❛r❣✐♥❛❧ ♣♦st❡r✐♦r ♠❡❛♥ ♦❢ ❜i ❛❢t❡r ✐♥t❡❣r❛t✐♥❣ ❛❝r♦ss t❤❡ ♣♦st❡r✐♦r ♦❢ β ✉♥❞❡r ❛ ✈❛❣✉❡ ♣r✐♦r ❢♦r β ✭s❡❡ ❲❛❦❡✜❡❧❞ ♣❛❣❡ ✸✽✵✮
✶✶✵ ✴ ✶✻✵
i Pi❩ i●
✶✶✶ ✴ ✶✻✵
✵ )
✐✐❞
b ∼ ●❛(c, d)
e ∼ ●❛(c, d)
✶✶✷ ✴ ✶✻✵
ij β + bi + eij, i = ✶, . . . , n; j = ✶, . . . , ni
iid
b) ❛♥❞ eij iid
e)✳
i=✶ ni ♦❜s❡r✈❛t✐♦♥s✱
✶✶✸ ✴ ✶✻✵
✶✶✹ ✴ ✶✻✵
ni×ni
b + σ✷ e
b
b
b
b + σ✷ e
b
b
b
b
b + σ✷ e
b, σ✷ e)T ✐s t❤❡ ✈❡❝t♦r ♦❢ ✈❛r✐❛♥❝❡ ❝♦♠♣♦♥❡♥ts✳
b
b+σ✷ e ❛♠♦♥❣ ❛❧❧ ♣❛✐rs (Yij, Yik)✳
✶✶✺ ✴ ✶✻✵
✵ )✱ ✇❤❡r❡ β✵ ❛♥❞ ❚ ✵ ❛r❡ ❦♥♦✇♥ q✉❛♥t✐t✐❡s
e ✐s ■● ✇✐t❤ ❦♥♦✇♥ s❤❛♣❡ ❛♥❞ s❝❛❧❡ ♣❛r❛♠❡t❡rs c ❛♥❞ d✱
e; c, d) d
e)−(c+✶) ❡①♣
e
e ∼ ●❛(c, d) ✇❤❡r❡ d ✐s ❛ r❛t❡
b ∼ ●❛(g, h) ❢♦r ❦♥♦✇♥ s❤❛♣❡ ❛♥❞ r❛t❡
✶✶✻ ✴ ✶✻✵
✶ β|②, r❡st ∼ ◆p(♠, ❱ )✱ ✇❤❡r❡
❱
p×p =
−✶✱ ✇❤❡r❡ τe = ✶/σ✷
e
♠
p×✶ = ❱
vi = ✶/(τb + niτe) mi = viτe✶T
ni (② i − ❳ iβ) = viτe
ni
j=✶(yij − ①T ij β)
= (✶ − wi)✵ + wi(¯ yi − ¯ ①T
i β)✱ ✇❤❡r❡ wi =
niσ✷
b
niσ✷
b + σ✷ e ✸ τe|②, r❡st ∼ ●❛
✶✶✼ ✴ ✶✻✵
########### # Priors # ########### beta0<-rep(0,p) # Prior mean for beta T0<-diag(.01,p) # Prior precision for beta c<-d<-.001 # Gamma hyperpriors for taue, taub ######### # Inits # ######### taue<-1 # Error precision = 1/sigma2 b<-rep(0,n) # Random effects taub<-1 # Random Effects precision ################# # GIBBS SAMPLER # ################# for (i in 1:nsim) { # Update beta vbeta<-solve(T0+taue*crossprod(X,X)) mbeta<-vbeta%*%(T0%*%beta0 + taue*crossprod(X,y-rep(b,nis))) beta<-c(rmvnorm(1,mbeta,vbeta)) # Update b vb<-1/(taub+nis*taue) mb<-vb*(taue*tapply(y-X%*%beta,id,sum)) # tapply sums (y-xbeta)'s for each subject b<-rnorm(n,mb,sqrt(vb)) # Update taue tmp<-d+crossprod(y-X%*%beta-rep(b,nis))/2 taue<-rgamma(1,c+N/2,tmp) # Update taub tmp<-c(d+crossprod(b)/2) taub<-rgamma(1,c+n/2,tmp) ################# # Store Results # ################# if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Sigmae[j]<-1/taue Sigmab[j]<-1/taub } if (i%%100==0) print(i) # 2.6 seconds to run 1000 iterations with n=1000 and N=5525 }
✶✶✽ ✴ ✶✻✵
iid
b)
iid
e)
e
b
† ❋r♦♠ ❙❆❙ Pr♦❝ ◆▲▼■❳❊❉✳ ✶✶✾ ✴ ✶✻✵
✺❖✬❍❛❣❛♥ ❛♥❞ ▲❡♦♥❛r❞ ✭✶✾✼✺✮❀ ❆③③❛❧✐♥✐✱ ✶✾✽✺
✶✷✵ ✴ ✶✻✵
d
✶✷✶ ✴ ✶✻✵
−6 −4 −2 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
SN Densities
x f(x) N(0,1) SN(0,1,3) SN(0,1,−3) SN(0,2,−3) ✶✷✷ ✴ ✶✻✵
✶✷✸ ✴ ✶✻✵
ij β + bi + eij, i = ✶, . . . , n; j = ✶, . . . , ni,
iid
b) ✇✐t❤ σ✷ b = ω✷/(✶ + α✷)✳
✶✷✹ ✴ ✶✻✵
✵
b)
ψ )
b ∼ ●❛(g, h)
✶✷✺ ✴ ✶✻✵
❚❤❡ ●✐❜❜ s❛♠♣❧❡r ❢♦r t❤❡ ❙◆ r❛♥❞♦♠ ✐♥t❡r❝❡♣t ♠♦❞❡❧ ✐s✿
✶ ❉r❛✇ β ❢r♦♠ ◆p(♠, ❱ )✱ ✇❤❡r❡
❱
p×p =
−✶ ♠
p×✶ = ❱
v = ✶/(✶ + τbψ✷) mi = vτbψbi
✸ ❋♦r i = ✶, . . . , n✱ ✉♣❞❛t❡ ψ ❢r♦♠ ◆(m, v)✱ ✇❤❡r❡
v = ✶/(τψ + τb✇ T✇) ✇❤❡r❡ ✇ = (w✶, . . . , wn)T m = v(τψµψ + τb✇ T❜)
✹ ❋♦r i = ✶, . . . , n✱ ❞r❛✇ bi ❢r♦♠ ◆(mi, vi)✱ ✇❤❡r❡
vi = ✶/(τb + niτe) mi = vi
ni (② i − ❳ iβ)
✷ , d + (② − ❳β − ❩❜)T(② − ❳β − ❩❜) ✷
✷, h + (❜ − ψ✇)T(❜ − ψ✇)
α = ψ√τb ω =
❙❡❡ ❙◆ ❘❛♥❞♦♠ ■♥t❡r❝❡♣t✳r ❢♦r ❞❡t❛✐❧s ✶✷✻ ✴ ✶✻✵
for (i in 1:nsim){ # Update beta vbeta<-solve(T0+taue*crossprod(X,X)) mbeta<-vbeta%*%(T0%*%beta0 + taue*crossprod(X,y-rep(b,nis))) beta<-c(rmvnorm(1,mbeta,vbeta)) # Update w v<-1/(1+taub*psi^2) m<-v*taub*psi*b w<-rtnorm(n,m,sqrt(v),lower=0) # Update psi v<-1/(taupsi+taub*crossprod(w)) m<-v*(taupsi*mupsi+taub*crossprod(w,b)) psi<-rnorm(1,m,sqrt(v)) # Update b vb<-1/(taub+nis*taue) mb<-vb*(taub*psi*w+taue*tapply(y-X%*%beta,id,sum)) b<-rnorm(n,mb,sqrt(vb)) # Update taue taue<-rgamma(1,0.01+N/2,0.01+crossprod(y-X%*%beta-rep(b,nis))/2) # Update taub tmp<-c(0.01+crossprod(b-psi*w)/2) taub<-rgamma(1,0.01+n/2,tmp) ############################### # Transform and Store Results # ############################### if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Omega[j]<-sqrt(1/taub+psi^2) Alpha[j]<-psi*sqrt(taub) Sigmae[j]<-1/taue B[j,]<-b } if (i%%100==0) print(i) # 38 seconds to run 10K iterations }
✶✷✼ ✴ ✶✻✵
iid
iid
e)
e
† ❆ss✉♠✐♥❣ ♥♦r♠❛❧ r❛♥❞♦♠ ❡✛❡❝ts✳ ✶✷✽ ✴ ✶✻✵
200 400 600 800 1000 0.7 0.9 1.1 1.3 Iteration β1 200 400 600 800 1000 1.8 2.0 2.2 2.4 Iteration ω 200 400 600 800 1000 −3.5 −3.0 −2.5 −2.0 −1.5 Iteration α 200 400 600 800 1000 1.95 2.05 2.15 Iteration σe
✶✷✾ ✴ ✶✻✵
✶✸✵ ✴ ✶✻✵
4 6 8 10 10 20 30 40 50 time Y
Subject−Specific Trajectories Indivdual Data Points
✶✸✶ ✴ ✶✻✵
ni×✶ = ❳ i ni×p β p×✶ + ❩ i ni×q ❜i q×✶ + ❡i ni×✶,
✶✸✷ ✴ ✶✻✵
N×✶ = ❳ N×p β p×✶ + ❩ N×qn ❜ qn×✶ + ❡ N×✶,
✶✸✸ ✴ ✶✻✵
N×✷n =
✶✸✹ ✴ ✶✻✵
ij β + b✶i + tijb✷i + eij,
✶
✷
ij β
✶ + ✷tijσ✶✷ + t✷ ijσ✷ ✷ + σ✷ e
✶ + (tij + tik)σ✶✷ + tijtikσ✷ ✷
✻◆♦t❡ ❝❤❛♥❣❡ ♦❢ ♥♦t❛t✐♦♥ ❢r♦♠ ● t♦ Σb
✶✸✺ ✴ ✶✻✵
✶✸✻ ✴ ✶✻✵
ni×✶
ni×ni
ni×q Σb q×q ❩ T i q×ni
e■ ni
i + ❘i = Σi, s❛②✳
✶✸✼ ✴ ✶✻✵
✵ )
b
✵ )
✶✸✽ ✴ ✶✻✵
✷tr(❈ ✵Σ−✶),
✶✸✾ ✴ ✶✻✵
iid
i=✶ ③i③T i ✳ ❚❤❡♥✱
✵
✵),
✶✹✵ ✴ ✶✻✵
❖♥❡ ❝❛♥ s❤♦✇ t❤❛t t❤❡ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧s ❢♦r t❤❡ ❧✐♥❡❛r r❛♥❞♦♠ s❧♦♣❡ ♠♦❞❡❧ ❛r❡✿
✶ β|②, r❡st ∼ ◆p(♠, ❱ )✱ ✇❤❡r❡
❱
p×p =
−✶ ♠
p×✶ = ❱
❱
✷n×✷n = (■ n ⊗ Σ−✶ b
+ τe❩ T❩)−✶ ♠
✷n×✶ = τe❱ ❩ T(② − ❳β) ✸ τe|②, r❡st ∼ ●❛
✷ , d + (② − ❳β − ❩❜)T(② − ❳β − ❩❜) ✷
n×✷ =
b✶✶ b✷✶ ✳ ✳ ✳ ✳ ✳ ✳ b✶n b✷n ✭❡①t❡♥s✐♦♥s t♦
q > ✷ ❛r❡ str❛✐❣❤t❢♦r✇❛r❞✮
❙❡❡ ❘❛♥❞♦♠ ❙❧♦♣❡✳r ❢♦r ❞❡t❛✐❧s
✶✹✶ ✴ ✶✻✵
########### # Priors # ########### beta0<-rep(0,p) # Prior Mean for beta T0<-diag(.01,p) # Prior Precision Matrix of beta (vague), independent d0<-g0<-.001 # Hyperpriors for tau nu0<-3 # DF for Wishart prior on Sigmab (G in notes) C0<-diag(2) # Scale matrix for IW Prior on Sigmab d<-d0+N/2 # Posterior df for taue nu<-nu0+n # Posterior df for Sigmab ################# # GIBBS SAMPLER # ################# for (i in 1:nsim) { # Update Beta vbeta<-solve(T0+taue*crossprod(X,X)) mbeta<-vbeta%*%(T0%*%beta0 + taue*crossprod(X,y-Z%*%b)) beta<-c(rmvnorm(1,mbeta,vbeta)) # Update b precb<-diagn%x%Taub+taue*crossprod(Z) # Posterior precision mb<-taue*crossprod(Z,y-X%*%beta) # Likelihood contribution to posterior mean b<-rmvnorm.canonical(1,mb,precb)[1,] # Update without inverting using spam package btmp[,]<-b # More efficient to fill in a pre-defined matrix object bmat<-t(btmp) # Update taue zb<-Z%*%b g<-g0+crossprod(y-X%*%beta-zb)/2 taue<-rgamma(1,d,as.numeric(g)) # Update Taub Sigma.b<-riwish(nu,C0+crossprod(bmat,bmat)) Taub<-solve(Sigma.b) ################# # Store Results # ################# if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Sigmae[j]<-1/taue Sigmab[j,]<-c(Sigma.b) } if (i%%50==0) print(i) # Very slow! } 1
✶✹✷ ✴ ✶✻✵
✶✹✸ ✴ ✶✻✵
✶
✷
✶
✷
✶|✷) = ◆
✶(✶ − ρ✷)
✷|✶) = ◆
✷(✶ − ρ✷)
✶(✶ − ρ✷)■ n
✷(✶ − ρ✷)■ n
✶✹✺ ✴ ✶✻✵
❚❤✐s ❧❡❛❞s t♦ ❛ ♠♦❞✐✜❡❞ ●✐❜❜s s❛♠♣❧❡r✿
✶ ❯♣❞❛t❡ β✱ τe✱ Σb ❛s ❜❡❢♦r❡ ✷ ❋r♦♠ Σb✱ ❢♦r♠
ρ = σ✶✷ σ✶σ✷ τ✶|✷ = ✶ σ✷
✶(✶ − ρ✷)
τ✷|✶ = ✶ σ✷
✷(✶ − ρ✷) ✸ ❋♦r i = ✶, . . . , n✱ ✉♣❞❛t❡ b✶i ❣✐✈❡♥ b✷i ❢r♦♠ ◆(mi, vi)✱ ✇❤❡r❡
vi = ✶ τ✶|✷ + τeni ✭♥♦t❡ t❤❡ s✐♠✐❧❛r✐t② t♦ t❤❡ r❛♥❞♦♠ ✐♥t❡r❝❡♣t ♠♦❞❡❧✮ mi = vi[τ✶|✷µ✶|✷ + τe
ni
(yij − ①T
ij β − tijb✷i)]
µ✶|✷ = ρσ✶ σ✷ b✷i
✹ ❋♦r i = ✶, . . . , n✱ ✉♣❞❛t❡ b✷i ❣✐✈❡♥ b✶i ❢r♦♠ ◆(mi, vi)✱ ✇❤❡r❡
vi = ✶ τ✷|✶ + τe ni
j=✶ t✷ ij
mi = vi[τ✷|✶µ✷|✶ + τe
ni
tij(yij − ①T
ij β − b✶i)]
µ✷|✶ = ρσ✷ σ✶ b✶i ❙❡❡ ❘❛♥❞♦♠ ❙❧♦♣❡ ❈♦♥❞✐t✐♦♥❛❧✳r ❢♦r ❞❡t❛✐❧s ✶✹✻ ✴ ✶✻✵
# Inits taue<-tau12<-tau21<-1 # Error precision = 1/sigma2, conditional precision of b1|b2, etc sigmab1<-sigmab2<-1 # Marginal variances of b1 and b2 rhob<-0 # Corr(b1,b2) b1<-b2<-rep(0,n) # Random effects (int and slope) beta<-rep(0,p) # Posterior mean and var of beta # Gibbs for (i in 1:nsim) { # Update beta vbeta<-solve(prec0+taue*crossprod(X,X)) mbeta<-vbeta%*%(prec0%*%beta0 + taue*crossprod(X,y-rep(b1,nis)-t*rep(b2,nis))) beta<-c(rmvnorm(1,mbeta,vbeta)) # Update b1|b2 vb<-1/(tau12+taue*nis) mu12<-rhob*sqrt(sigmab1/sigmab2)*b2 # Prior Mean of b1|b2 mb<-vb*(tau12*mu12+taue*tapply(y-X%*%beta-t*rep(b2,nis),id,sum)) b1<-rnorm(n,mb,sqrt(vb)) # Update b2|b1 vb<-1/(tau21+taue*tapply(t^2,id,sum)) mu21<-rhob*sqrt(sigmab2/sigmab1)*b1 # Prior Mean of b2|b1 mb<-vb*(tau21*mu21+taue*tapply(t*(y-X%*%beta-rep(b1,nis)),id,sum)) b2<-rnorm(n,mb,sqrt(vb)) # Update taue g<-g0+crossprod(y-X%*%beta-rep(b1,nis)-t*rep(b2,nis))/2 taue<-rgamma(1,d0+N/2,g) # Update Sigma.b bmat<-cbind(b1,b2) Sigma.b<-riwish(nu0+n,C0+crossprod(bmat)) sigmab1<-Sigma.b[1,1] # Marginal variance of b1 sigmab2<-Sigma.b[2,2] # Marginal variance of b2 rhob<-Sigma.b[1,2]/sqrt(sigmab1*sigmab2) # Corr(b1,b2) tau12<-1/(sigmab1*(1-rhob^2)) # Conditional precision of b1|b2 tau21<-1/(sigmab2*(1-rhob^2)) # Conditional precision of b2|b1 ################# # Store Results # ################# if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Sigmae[j]<-1/taue Sigmab[j,]<-c(Sigma.b) B1[j,]<-b1[1:10] # Store first 10 random effects B2[j,]<-b2[1:10] } } # 66 seconds to run 10,000 iterations with n=1000 1
✶✹✼ ✴ ✶✻✵
iid
iid
e)
❚❛❜❧❡ ✶✶✿ ❘❡s✉❧ts ❢r♦♠ ❈♦♥❞✐t✐♦♥❛❧ ●✐❜❜s ❆❧❣♦r✐t❤♠✳
✶
✷
e
† ❋r♦♠ ❙❆❙ Pr♦❝ ▼✐①❡❞✳
✶✹✽ ✴ ✶✻✵
✶✹✾ ✴ ✶✻✵
ij β + bi
b
✶ ❋♦r ❛❧❧ (i, j✮✱ ✉♣❞❛t❡ t❤❡ P● ✇❡✐❣❤ts✱ ωij✱ t♦ ✐♥❝❧✉❞❡ bi ✷ ❋♦r♠ t❤❡ ❝♦rr❡s♣♦♥❞✐♥❣ ❧❛t❡♥t ♥♦r♠❛❧ ✈❛r✐❛❜❧❡s✱ Zij ✸ ❈♦♥❞✐t✐♦♥❛❧ ♦♥ t❤❡ ❧❛t❡♥t ♥♦r♠❛❧ ✈❛r✐❛❜❧❡s✱ ✉♣❞❛t❡ β ❛♥❞ bi
✶✺✵ ✴ ✶✻✵
❙♣❡❝✐✜❝❛❧❧②✱ t❤❡ ●✐❜❜s s❛♠♣❧❡r ♣r♦❝❡❡❞s ❛s ❢♦❧❧♦✇s✿
✶ ❋♦r ❛❧❧ (i, j)✱ ✉♣❞❛t❡ ωij ❢r♦♠ P●(✶, ηij)✱ ✇❤❡r❡ ηij = ①T ij β + bi ✷ ❋♦r ❛❧❧ (i, j)✱ ❞❡✜♥❡ zij = yij−✶/✷ ωij ✸ ❯♣❞❛t❡ β ❢r♦♠ ◆p(♠, ❱ )✱ ✇❤❡r❡
❱ =
−✶ ♠ = ❱ ❚ ✵β✵ + ❳ T
p×N Ω N×N (③ − ❲ ❜)
, ✇❤❡r❡ Ω
N×N = ❞✐❛❣(ωij) ❛♥❞ ❲ N×n ❞❡♥♦t❡s t❤❡ r❛♥❞♦♠ ❡✛❡❝ts ❞❡s✐❣♥
♠❛tr✐①✳
✹ ❋♦r ❛❧❧ i✱ ✉♣❞❛t❡ bi ❢r♦♠ ◆(mi, vi)✱ ✇❤❡r❡
vi = ✶ τb + ni
j=✶ ωij
mi = vi ❲ T
i ✶×ni
Ωi
ni×ni
(③i − ❳ iβ)
= vi
ni
ωij(zij − ①T
ij β) ✺ ❯♣❞❛t❡ τb ❢r♦♠ ●❛
✷, d + ❜T❜/✷
❤②♣❡r♣r✐♦rs ❙❡❡ ▲♦❣✐st✐❝ ❘❛♥❞♦♠ ■♥t❡r❝❡♣t✳r ❢♦r ❞❡t❛✐❧s ✶✺✶ ✴ ✶✻✵
# Priors beta0<-rep(0,p) # Prior mean for beta T0<-diag(.01,p) # Prior precision for beta c<-d<-0.01 # Gamma hyperpriors # Inits beta<-rep(0,p) b<-rep(0,n) # Random effect taub<-1 # Random effect precision ######### # Gibbs # ######### tmp<-proc.time() for (i in 1:nsim){ # Update z mu<-X%*%beta+rep(b,nis)
z<-(y-1/2)/omega # Update beta v<-solve(crossprod(sqrt(omega)*X)+T0) m<-v%*%(T0%*%beta0+t(sqrt(omega)*X)%*%(sqrt(omega)*(z-rep(b,nis)))) beta<-c(rmvnorm(1,m,v)) # Update b vb<-1/(taub+tapply(omega,id,sum)) mb<-vb*(tapply(omega*(z-X%*%beta),id,sum)) b<-rnorm(n,mb,sqrt(vb)) # Update taub taub<-rgamma(1,c+n/2,d+crossprod(b)/2) # Store if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Sigmab[j]<-1/taub } if (i%%100==0) print(i) } tot.time<-proc.time()-tmp # 8.5 secs to run 1000 iterations
✶✺✷ ✴ ✶✻✵
✶ ❋♦r ❛❧❧ (i, j)✱ ❞r❛✇ ωij ❢r♦♠ ✐ts P●(yij + r, ηi) ❞✐str✐❜✉t✐♦♥✱
i β + bi ✷ ❋♦r ❛❧❧ (i, j)✱ ❞❡✜♥❡ zij = yij − r
✸ ❯♣❞❛t❡ β ❢r♦♠ ✐ts ◆p(♠, ❱ ) ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧✱ ✇❤❡r❡
✹ ❯♣❞❛t❡ ❜ ❛♥❞ τb ❢r♦♠ t❤❡✐r ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧s ❛♥❛❧♦❣♦✉s t♦
✺ ❯♣❞❛t❡ r ✉s✐♥❣ r❛♥❞♦♠✲✇❛❧❦ ▼❍ ♦r ❛ ❝♦♥❥✉❣❛t❡ ●❛♠♠❛
✶✺✸ ✴ ✶✻✵
# Priors similar to fixed effects NB model, but with gamma for taub # Inits beta<-rep(0,p) b<-rep(0,n) # Random effects taub<-1 # Random effect precision r<-1 # Inverse Dispersion Acc<-0 # Acceptance rate indicator s<-.005 # Proposal variance ######## # MCMC # ######## for (i in 1:nsim){ # Update z and beta eta<-X%*%beta+rep(b,nis)
# Polya weights z<-(y-r)/(2*omega) v<-solve(crossprod(X*sqrt(omega))+T0) m<-v%*%(T0%*%beta0+t(sqrt(omega)*X)%*%(sqrt(omega)*(z-rep(b,nis)))) beta<-c(rmvnorm(1,m,v)) # Update b vb<-1/(taub+tapply(omega,id,sum)) mb<-vb*(tapply(omega*(z-X%*%beta),id,sum)) b<-rnorm(n,mb,sqrt(vb)) # Update taub taub<-rgamma(1,c+n/2,d+crossprod(b)/2) # Update r q<-1/(1+exp(eta)) # dnegbinom uses q=1-p rnew<-rtnorm(1,r,sqrt(s),lower=0) ratio<-sum(dnbinom(y,rnew,q,log=T))-sum(dnbinom(y,r,q,log=T))+ dtnorm(rnew,r,sqrt(s),0,log=T)-dtnorm(r,rnew,sqrt(s),0,log=T) if (log(runif(1))<ratio) { r<-rnew Acc<-Acc+1 } if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Sigmab[j]<-1/taub R[j]<-r } if (i%%100==0) print(i) # 23 secs to run 1000 iterations with n=1000 and N=5515 }
✶✺✹ ✴ ✶✻✵
# Priors similar to fixed effects NB gibbs model, but with gamma for taub # Inits beta<-rep(0,p) b<-rep(0,n) # Random effects taub<-1 # Random effect precision r<-1 # Inverse Dispersion ######## # MCMC # ######## for (i in 1:nsim){ # Update z eta<-X%*%beta+rep(b,nis)
# Polya weights z<-(y-r)/(2*omega) # Update beta v<-solve(crossprod(X*sqrt(omega))+T0) m<-v%*%(T0%*%beta0+t(sqrt(omega)*X)%*%(sqrt(omega)*(z-rep(b,nis)))) beta<-c(rmvnorm(1,m,v)) # Update b vb<-1/(taub+tapply(omega,id,sum)) mb<-vb*(tapply(omega*(z-X%*%beta),id,sum)) b<-rnorm(n,mb,sqrt(vb)) # Update taub taub<-rgamma(1,c+n/2,d+crossprod(b)/2) # 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 psi eta<-X%*%beta+rep(b,nis) psi<-exp(eta)/(1+exp(eta)) r<-rgamma(1,a+sum(l),b-sum(log(1-psi))) if (i> burn & i%%thin==0) { j<-(i-burn)/thin Beta[j,]<-beta Sigmab[j]<-1/taub R[j]<-r B[j,]<-b } if (i%%100==0) print(i) # 49 secs to run 1000 iterations with n=1000 and N=5515
✶✺✺ ✴ ✶✻✵
ij β + b✶i
ij γ + b✷i + eij
✶
✷
e ) i = ✶, . . . , n; j = ✶, . . . , ni
✶✺✻ ✴ ✶✻✵
✶ ❆✉❣♠❡♥t t❤❡ ❆❧❜❡rt ✫ ❈❤✐❜ s❛♠♣❧❡r ✇✐t❤ ❛ r❛♥❞♦♠
✷ ❋✐t ❛ ❧✐♥❡❛r r❛♥❞♦♠ ✐♥t❡r❝❡♣t ♠♦❞❡❧ t♦ ✉♣❞❛t❡ t❤❡ ❧✐♥❡❛r
✸ ❯♣❞❛t❡ t❤❡ r❛♥❞♦♠ ✐♥t❡r❝❡♣ts ✉s✐♥❣ t❤❡ ❝♦♥❞✐t✐♦♥❛❧ ❢♦r♠✉❧❛
✹ ❯♣❞❛t❡ Σb ❢r♦♠ ✐ts ■❲ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧
✶✺✼ ✴ ✶✻✵
✶ Pr♦❜✐t ♠♦❞❡❧✿
ij β + bi, ✶)
✷ ▲✐♥❡❛r ♠♦❞❡❧✿
✸ ❯♣❞❛t❡ Σb ❢r♦♠ ✐ts ❝♦♥❥✉❣❛t❡ ■❲ ❢✉❧❧ ❝♦♥❞✐t✐♦♥❛❧ ❛♥❞ r❡tr✐❡✈❡ τ✶|✷✱
✶✺✽ ✴ ✶✻✵
# Priors same as for probit and linear model with IW prior for (b1,b2) vbeta<-solve(T0+crossprod(X,X)) # Posterior Var(beta) -- can do outside loop # GIBBS SAMPLER for (i in 1:nsim) { # Binomial Model # Draw Latent Variable, z muz<-X%*%beta+rep(b1,nis) # Mean of z z[y1==0]<-qnorm(runif(N,0,pnorm(0,muz)),muz)[y1==0] z[y1==1]<-qnorm(runif(N,pnorm(0,muz),1),muz)[y1==1] # Update beta mbeta <- vbeta%*%(T0%*%beta0+crossprod(X,z-rep(b1,nis))) beta<-c(rmvnorm(1,mbeta,vbeta)) # Update b1|b2 taub12<-taub1/(1-rhob^2) # Prior precision of b1|b2 mb12<-rhob*sqrt(taub2/taub1)*b2 # Prior mean of b1|b2 vb1<-1/(nis+taub12) # Posterior var of b1|b2,y mb1<-vb1*(taub12*mb12+tapply(z-X%*%beta,id,sum)) # Posterior mean of b1|b2,y b1<-rnorm(n,mb1,sqrt(vb1)) # Linear Model # Update gamma vgam<-solve(G0+taue*crossprod(X)) # G0 is prior precision of gamma mgam<-vgam%*%(G0%*%gamma0 + taue*crossprod(X,y2-rep(b2,nis))) gamma<-c(rmvnorm(1,mgam,vgam)) # Update taue l<-l0+crossprod(y2-X%*%gamma-rep(b2,nis))/2 taue<-rgamma(1,d0+N/2,l) # Update b2|b1 taub21<-taub2/(1-rhob^2) # Prior precision of b2|b1 mb21<-rhob*sqrt(taub1/taub2)*b1 # Prior mean of b2|b1 vb2<-1/(taue*nis+taub21) # Posterior var of b2|b1 mb2<-vb2*(taub21*mb21+taue*tapply(y2-X%*%gamma,id,sum)) b2<-rnorm(n,mb2,sqrt(vb2)) b<-cbind(b1,b2) # Update variance components Sigmab<-riwish(nu0+n,c0+crossprod(b)) rhob<-Sigmab[1,2]/sqrt(Sigmab[1,1]*Sigmab[2,2]) taub1<-1/Sigmab[1,1] taub2<-1/Sigmab[2,2] } # 7.4 seconds to run 1000 iterations with n=1000 and N=5441
✶✺✾ ✴ ✶✻✵
✶✻✵ ✴ ✶✻✵