Performing Bayesian analysis in Stata using WinBUGS Tom Palmer, John - - PowerPoint PPT Presentation

performing bayesian analysis in stata using winbugs
SMART_READER_LITE
LIVE PREVIEW

Performing Bayesian analysis in Stata using WinBUGS Tom Palmer, John - - PowerPoint PPT Presentation

Performing Bayesian analysis in Stata using WinBUGS Tom Palmer, John Thompson & Santiago Moreno Department of Health Sciences, University of Leicester, UK 13 th UK Stata Users Group Meeting, 10 September 2007 Tom Palmer (Leicester) Running


slide-1
SLIDE 1

Performing Bayesian analysis in Stata using WinBUGS

Tom Palmer, John Thompson & Santiago Moreno

Department of Health Sciences, University of Leicester, UK

13th UK Stata Users Group Meeting, 10 September 2007

Tom Palmer (Leicester) Running WinBUGS from Stata 1 / 24

slide-2
SLIDE 2

Outline

1

The Bayesian approach & WinBUGS

2

The winbugsfromstata package

3

How to run an analysis

4

Summary & developments

Tom Palmer (Leicester) Running WinBUGS from Stata 2 / 24

slide-3
SLIDE 3

The Bayesian approach

Bayes Theorem

Posterior ∝ Likelihood × prior

Tom Palmer (Leicester) Running WinBUGS from Stata 3 / 24

slide-4
SLIDE 4

The Bayesian approach

Bayes Theorem

Posterior ∝ Likelihood × prior Direct probability statements - not frequentist - subjective

Tom Palmer (Leicester) Running WinBUGS from Stata 3 / 24

slide-5
SLIDE 5

The Bayesian approach

Bayes Theorem

Posterior ∝ Likelihood × prior Direct probability statements - not frequentist - subjective Complex posterior marginal distributions - estimation via simulation

Tom Palmer (Leicester) Running WinBUGS from Stata 3 / 24

slide-6
SLIDE 6

The Bayesian approach

Bayes Theorem

Posterior ∝ Likelihood × prior Direct probability statements - not frequentist - subjective Complex posterior marginal distributions - estimation via simulation Markov chain Monte Carlo (MCMC) methods

Tom Palmer (Leicester) Running WinBUGS from Stata 3 / 24

slide-7
SLIDE 7

WinBUGS

Bayesian statistics using Gibbs sampling

Tom Palmer (Leicester) Running WinBUGS from Stata 4 / 24

slide-8
SLIDE 8

WinBUGS

Bayesian statistics using Gibbs sampling MRC Biostatistics unit

http://www.mrc-bsu.cam.ac.uk/bugs

Tom Palmer (Leicester) Running WinBUGS from Stata 4 / 24

slide-9
SLIDE 9

WinBUGS

Bayesian statistics using Gibbs sampling MRC Biostatistics unit

http://www.mrc-bsu.cam.ac.uk/bugs

Health Economics, Medical Statistics

Tom Palmer (Leicester) Running WinBUGS from Stata 4 / 24

slide-10
SLIDE 10

WinBUGS

Bayesian statistics using Gibbs sampling MRC Biostatistics unit

http://www.mrc-bsu.cam.ac.uk/bugs

Health Economics, Medical Statistics Disadvantages: data management, post-processing of results, graphics

Tom Palmer (Leicester) Running WinBUGS from Stata 4 / 24

slide-11
SLIDE 11

The winbugsfromstata package

Stata interface to WinBUGS [Thompson et al., 2006]

http://www2.le.ac.uk/departments/health-sciences/extranet/ BGE/genetic-epidemiology/gedownload/information

Tom Palmer (Leicester) Running WinBUGS from Stata 5 / 24

slide-12
SLIDE 12

The winbugsfromstata package

Tom Palmer (Leicester) Running WinBUGS from Stata 6 / 24

slide-13
SLIDE 13

How to run an analysis

Tom Palmer (Leicester) Running WinBUGS from Stata 7 / 24

slide-14
SLIDE 14

help winbugs

Tom Palmer (Leicester) Running WinBUGS from Stata 8 / 24

slide-15
SLIDE 15

Example analysis: Schools

Schools example [Goldstein et al., 1993],[Spiegelhalter et al., 2004]

Tom Palmer (Leicester) Running WinBUGS from Stata 9 / 24

slide-16
SLIDE 16

Example analysis: Schools

Schools example [Goldstein et al., 1993],[Spiegelhalter et al., 2004] Between-school variation in exam results from inner London schools

Tom Palmer (Leicester) Running WinBUGS from Stata 9 / 24

slide-17
SLIDE 17

Example analysis: Schools

Schools example [Goldstein et al., 1993],[Spiegelhalter et al., 2004] Between-school variation in exam results from inner London schools Standardized mean scores (Y ) 1,978 pupils, 38 schools

Tom Palmer (Leicester) Running WinBUGS from Stata 9 / 24

slide-18
SLIDE 18

Example analysis: Schools

Schools example [Goldstein et al., 1993],[Spiegelhalter et al., 2004] Between-school variation in exam results from inner London schools Standardized mean scores (Y ) 1,978 pupils, 38 schools LRT: London Reading Test, VR: verbal reasoning, Gender intake of school, denomination of school

Tom Palmer (Leicester) Running WinBUGS from Stata 9 / 24

slide-19
SLIDE 19

Data for the Schools example

−2 2 4 −2 2 4 −2 2 4 −2 2 4 −2 2 4 −2 2 4 −40 −20 20 40 −40 −20 20 40 −40 −20 20 40 −40 −20 20 40 −40 −20 20 40 −40 −20 20 40 −40 −20 20 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

girl boy Y LRT

Graphs by school Tom Palmer (Leicester) Running WinBUGS from Stata 10 / 24

slide-20
SLIDE 20

The model

Hierarchical model; specified the mean and variance

Tom Palmer (Leicester) Running WinBUGS from Stata 11 / 24

slide-21
SLIDE 21

The model

Hierarchical model; specified the mean and variance Model: Yij ∼ N(µij, τij) µij = γ1j + γ2jLRT ij + γ3jVR1ij + β1LRT 2

ij + β2VR2ij

+ β3Girlij + β4Gschj + β5Bschj + β6CEschj + β7RCschj + β8Oschj log τij = θ + φLRT ij

Tom Palmer (Leicester) Running WinBUGS from Stata 11 / 24

slide-22
SLIDE 22

WinBUGS model statement

model{ for(p in 1 : N){ Y[p] ~ dnorm(mu[p], tau[p]) mu[p] <- alpha[school[p], 1] + alpha[school[p], 2] * LRT[p] + alpha[school[p], 3] * VR[p, 1] + beta[1] * LRT2[p] + beta[2] * VR[p, 2] + beta[3] * Gender[p] + beta[4] * School.gender[p, 1] + beta[5] * School.gender[p, 2] + beta[6] * School.denom[p, 1] + beta[7] * School.denom[p, 2] + beta[8] * School.denom[p, 3] log(tau[p]) <- theta + phi * LRT[p] sigma2[p] <- 1 / tau[p] LRT2[p] <- LRT[p] * LRT[p] } min.var <- exp(-(theta + phi * (-34.6193))) # lowest LRT score = -34.6193 max.var <- exp(-(theta + phi * (37.3807))) # highest LRT score = 37.3807 # Priors for fixed effects: for (k in 1 : 8){ beta[k] ~ dnorm(0.0, 0.0001) } theta ~ dnorm(0.0, 0.0001) phi ~ dnorm(0.0, 0.0001) # Priors for random coefficients: for (j in 1 : M) { alpha[j, 1 : 3] ~ dmnorm(gamma[1:3 ], T[1:3 ,1:3 ]) alpha1[j] <- alpha[j,1] } # Hyper-priors: gamma[1 : 3] ~ dmnorm(mn[1:3 ], prec[1:3 ,1:3 ]) T[1 : 3, 1 : 3 ] ~ dwish(R[1:3 ,1:3 ], 3) }

Tom Palmer (Leicester) Running WinBUGS from Stata 12 / 24

slide-23
SLIDE 23

Do-file for the example

// winbugsfromstata demo, 16august2007 cd "Z:/conferences/stata.users.uk.2007/schools" wbdecode, file(Schoolsdata.txt) clear wbscript, sav(‘c(pwd)’/script.txt, replace) /// model(‘c(pwd)’/Schoolsmodel.txt) /// data(‘c(pwd)’/Schoolsdata.txt) /// inits(‘c(pwd)’/Schoolsinits.txt) /// coda(‘c(pwd)’/out) /// burn(500) update(1000) /// set(beta gamma phi theta) dic /// log(‘c(pwd)’/winbugslog.txt) /// quit wbrun , sc(‘c(pwd)’/script.txt) /// win(Z:/winbugs/WinBUGS14/WinBUGS14.exe) clear set memory 500m wbcoda, root(out) clear wbstats gamma* beta* phi theta wbtrace beta_1 gamma_1 phi theta wbdensity beta_1 gamma_1 phi theta wbac beta_1 gamma_1 phi theta wbhull beta_1 beta_2 gamma_2, peels(1 5 10 25) wbgeweke beta_1 gamma_1 phi theta wbdic using winbugslog.txt

Tom Palmer (Leicester) Running WinBUGS from Stata 13 / 24

slide-24
SLIDE 24

wbrun screenshot 1

Tom Palmer (Leicester) Running WinBUGS from Stata 14 / 24

slide-25
SLIDE 25

wbrun screenshot 2

Tom Palmer (Leicester) Running WinBUGS from Stata 15 / 24

slide-26
SLIDE 26

Stata output

wbstats output

. wbstats gamma* beta* phi theta Parameter n mean sd sem median 95% CrI gamma_1 500

  • 0.715

0.103 0.0179

  • 0.715 (
  • 0.951,
  • 0.523 )

gamma_2 500 0.031 0.010 0.0005 0.031 ( 0.010, 0.052 ) gamma_3 500 0.967 0.105 0.0225 0.972 ( 0.750, 1.168 ) beta_1 500 0.000 0.000 0.0000 0.000 ( 0.000, 0.000 ) beta_2 500 0.433 0.072 0.0099 0.435 ( 0.284, 0.576 ) beta_3 500 0.173 0.048 0.0031 0.172 ( 0.085, 0.271 ) beta_4 500 0.151 0.141 0.0230 0.164 (

  • 0.156,

0.392 ) beta_5 500 0.091 0.105 0.0150 0.087 (

  • 0.094,

0.318 ) beta_6 500

  • 0.279

0.183 0.0279

  • 0.290 (
  • 0.618,

0.108 ) beta_7 500 0.170 0.105 0.0158 0.169 (

  • 0.029,

0.380 ) beta_8 500

  • 0.109

0.209 0.0376

  • 0.124 (
  • 0.485,

0.357 ) phi 500

  • 0.003

0.003 0.0002

  • 0.003 (
  • 0.009,

0.003 ) theta 500 0.579 0.032 0.0016 0.579 ( 0.513, 0.649 )

Tom Palmer (Leicester) Running WinBUGS from Stata 16 / 24

slide-27
SLIDE 27

Stata output

wbstats output

. wbstats gamma* beta* phi theta Parameter n mean sd sem median 95% CrI gamma_1 500

  • 0.715

0.103 0.0179

  • 0.715 (
  • 0.951,
  • 0.523 )

gamma_2 500 0.031 0.010 0.0005 0.031 ( 0.010, 0.052 ) gamma_3 500 0.967 0.105 0.0225 0.972 ( 0.750, 1.168 ) beta_1 500 0.000 0.000 0.0000 0.000 ( 0.000, 0.000 ) beta_2 500 0.433 0.072 0.0099 0.435 ( 0.284, 0.576 ) beta_3 500 0.173 0.048 0.0031 0.172 ( 0.085, 0.271 ) beta_4 500 0.151 0.141 0.0230 0.164 (

  • 0.156,

0.392 ) beta_5 500 0.091 0.105 0.0150 0.087 (

  • 0.094,

0.318 ) beta_6 500

  • 0.279

0.183 0.0279

  • 0.290 (
  • 0.618,

0.108 ) beta_7 500 0.170 0.105 0.0158 0.169 (

  • 0.029,

0.380 ) beta_8 500

  • 0.109

0.209 0.0376

  • 0.124 (
  • 0.485,

0.357 ) phi 500

  • 0.003

0.003 0.0002

  • 0.003 (
  • 0.009,

0.003 ) theta 500 0.579 0.032 0.0016 0.579 ( 0.513, 0.649 )

regress γ2: 0.030, 95% C.I. (0.026, 0.034)

Tom Palmer (Leicester) Running WinBUGS from Stata 16 / 24

slide-28
SLIDE 28

Stata output

wbgeweke output

. wbgeweke beta_1 Parameter: beta_1 first 10.0% (n=50) vs last 50.0% (n=250) Means (se) 0.0003 ( 0.0000) 0.0003 ( 0.0000) Autocorrelations 0.3736 0.4114 Mean Difference (se) 0.0000 ( 0.0000) z = 1.030 p = 0.3031

Tom Palmer (Leicester) Running WinBUGS from Stata 17 / 24

slide-29
SLIDE 29

Stata output

wbgeweke output

. wbgeweke beta_1 Parameter: beta_1 first 10.0% (n=50) vs last 50.0% (n=250) Means (se) 0.0003 ( 0.0000) 0.0003 ( 0.0000) Autocorrelations 0.3736 0.4114 Mean Difference (se) 0.0000 ( 0.0000) z = 1.030 p = 0.3031

wbdic output

. wbdic using winbugslog.txt DIC statistics 1 DIC Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes Dbar Dhat pD DIC Y 4466.330 4393.470 72.861 4539.190 total 4466.330 4393.470 72.861 4539.190

Tom Palmer (Leicester) Running WinBUGS from Stata 17 / 24

slide-30
SLIDE 30

wbtrace output

.0002 .0004 .0006 beta_1 2000 4000 6000 8000 10000 Order −1 −.8 −.6 −.4 −.2 gamma_1 2000 4000 6000 8000 10000 Order −.015−.01−.005 .005 .01 phi 2000 4000 6000 8000 10000 Order .45 .5 .55 .6 .65 .7 theta 2000 4000 6000 8000 10000 Order

Tom Palmer (Leicester) Running WinBUGS from Stata 18 / 24

slide-31
SLIDE 31

wbdensity output

1000 2000 3000 4000 Density −.0002 .0002 .0004 .0006 Estimate

beta_1

1 2 3 4 Density −1.2 −1 −.8 −.6 −.4 −.2 Estimate

gamma_1

50 100 150 Density −.015 −.01 −.005 .005 .01 Estimate

phi

5 10 15 Density .45 .5 .55 .6 .65 .7 Estimate

theta

Tom Palmer (Leicester) Running WinBUGS from Stata 19 / 24

slide-32
SLIDE 32

wbac output

0.00 0.10 0.20 0.30 0.40 Autocorrelations of beta_1 10 20 30 40 Lag

Bartlett’s formula for MA(q) 95% confidence bands

−0.20 0.000.200.400.600.80 Autocorrelations of gamma_1 10 20 30 40 Lag

Bartlett’s formula for MA(q) 95% confidence bands

−0.05 0.00 0.05 0.10 Autocorrelations of phi 10 20 30 40 Lag

Bartlett’s formula for MA(q) 95% confidence bands

−0.02 0.00 0.02 0.04 0.06 Autocorrelations of theta 10 20 30 40 Lag

Bartlett’s formula for MA(q) 95% confidence bands

Tom Palmer (Leicester) Running WinBUGS from Stata 20 / 24

slide-33
SLIDE 33

wbhull output

+ .0002 .0004 .0006 beta_1 .2 .3 .4 .5 .6 beta_2 + .0002 .0004 .0006 beta_1 −.02 .02 .04 .06 .08 gamma_2 + .2 .3 .4 .5 .6 beta_2 −.02 .02 .04 .06 .08 gamma_2

Tom Palmer (Leicester) Running WinBUGS from Stata 21 / 24

slide-34
SLIDE 34

Summary

WinBUGS - easy & flexible

Tom Palmer (Leicester) Running WinBUGS from Stata 22 / 24

slide-35
SLIDE 35

Summary

WinBUGS - easy & flexible winbugsfromstata - data preparation, analysis of MCMC output, graphics

Tom Palmer (Leicester) Running WinBUGS from Stata 22 / 24

slide-36
SLIDE 36

Summary

WinBUGS - easy & flexible winbugsfromstata - data preparation, analysis of MCMC output, graphics Prior distributions - controversial

Tom Palmer (Leicester) Running WinBUGS from Stata 22 / 24

slide-37
SLIDE 37

Summary

WinBUGS - easy & flexible winbugsfromstata - data preparation, analysis of MCMC output, graphics Prior distributions - controversial Check complex Stata models - vague prior distributions

Tom Palmer (Leicester) Running WinBUGS from Stata 22 / 24

slide-38
SLIDE 38

Summary

WinBUGS - easy & flexible winbugsfromstata - data preparation, analysis of MCMC output, graphics Prior distributions - controversial Check complex Stata models - vague prior distributions Fit complex models not possible in Stata

Tom Palmer (Leicester) Running WinBUGS from Stata 22 / 24

slide-39
SLIDE 39

Developments

Bayesian residuals and model checking [Lu et al., 2007]

Tom Palmer (Leicester) Running WinBUGS from Stata 23 / 24

slide-40
SLIDE 40

Developments

Bayesian residuals and model checking [Lu et al., 2007] Automate WinBUGS model statement

Tom Palmer (Leicester) Running WinBUGS from Stata 23 / 24

slide-41
SLIDE 41

Developments

Bayesian residuals and model checking [Lu et al., 2007] Automate WinBUGS model statement Mac users: WinBUGS runs under Darwine

Tom Palmer (Leicester) Running WinBUGS from Stata 23 / 24

slide-42
SLIDE 42

Developments

Bayesian residuals and model checking [Lu et al., 2007] Automate WinBUGS model statement Mac users: WinBUGS runs under Darwine OpenBUGS (version 3.0.1), WinBUGS (version 1.4.2)

http://mathstat.helsinki.fi/openbugs/

Tom Palmer (Leicester) Running WinBUGS from Stata 23 / 24

slide-43
SLIDE 43

References

Goldstein, H., Rasbash, J., Yang, M., Woodhouse, G., Pan, H., Nuttall, D., and Thomas, S. (1993). A multilevel analysis of school examination results. Oxford Review of Education, 19(4):425–433. Lu, G., Ades, A. E., Sutton, A. J., Cooper, N. J., Briggs, A. H., and Caldwell, D. M. (2007). Meta-analysis of mixed treatment comparisons at multiple follow-up times. Statistics in Medicine. in press. Spiegelhalter, D. J., Thomas, A., Best, N., and Lunn, D. (2004). WinBUGS User Manual, version 1.4.1. MRC Biostatistics Unit, Cambridge, UK. Thompson, J., Palmer, T., and Moreno, S. (2006). Bayesian Analysis in Stata using WinBUGS. The Stata Journal, 6(4):530–549.

Acknowledgements MRC Capacity Building PhD Studentship in Genetic Epidemiology

Tom Palmer (Leicester) Running WinBUGS from Stata 24 / 24