【好书推荐】Bayesian Methods in Health Economics

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

http://www.crcpress.com/product/isbn/9781439895559

Features

• Provides an overview of Bayesian methods for cost-effectiveness analysis, and includes all necessary background on economics and Bayesian statistics
• Presents three detailed case studies of the cost-effectiveness analysis of health care interventions
• Includes several worked examples to guide through the process of health economic evaluation
• Contains extensive coverage of the practice of making Bayesian analysis integrating software such as JAGS and R, specifically for the application of health economic analysis
• Systematically describes the methodological issues related to the application of Bayesian inference and decision process in health economics
• Designed as a reference for students and practitioners working in the field of health economic evaluations and medical statistics
  

The book is linked to code with which to replicate the examples, and an associated R package (BCEA) can be used in real applications to produce systematic health economic evaluations of Bayesian models.

Summary

Health economics is concerned with the study of the cost-effectiveness of health care interventions. This book provides an overview of Bayesian methods for the analysis of health economic data. After an introduction to the basic economic concepts and methods of evaluation, it presents Bayesian statistics using accessible mathematics. The next chapters describe the theory and practice of cost-effectiveness analysis from a statistical viewpoint, and Bayesian computation, notably MCMC. The final chapter presents three detailed case studies covering cost-effectiveness analyses using individual data from clinical trials, evidence synthesis and hierarchical models and Markov models. The text uses WinBUGS and JAGS with datasets and code available online.

本帖隐藏的内容

Bayesian Methods in Health Economics.zip (2.75 MB, 需要: 55 个论坛币) 本附件包括：

• Bayesian Methods in Health Economics.pdf

2014-11-28 11:31:06 上传

需要: 55 个论坛币  [购买]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享1 收藏74 回帖

关键词：health econ Economics Economic Bayesian Methods 好书推荐 Health

Bayesian Methods in Health EconomicsCheck out my book (published with CRC - available also on Amazon, in ebook format too)

Table of contents Preface Get a promotional code - save 20% when ordering online from CRC website Read a sample for free BCEA - an R package to run Bayesian health economic evaluations (used throughout the book and specifically in the examples) Some discussion of the book in the blog can be found here, here, here and here A list of typos and corrections is discussed here R/JAGS code to run the examples (all files at once here) Instructions(probably best to read before downloading the files) Utilities Utils.R: script containing some utility functions, for examples to draw traceplots of MCMC chains, or estimating the parameters of suitable distributions to obtain given values for its mean and standard deviation [Needed to run most of the other scripts] Chapter 2 MCMC.R: script to run the Gibbs sampling simulations and check convergence, as in Figure 2.10 modelNormal.R: script to run the analysis of the Normal model (pages 69-73) modelNormal.txt: JAGS code for the Bayesian model (pages 69-73) phbirths.dta: Dataset used for the Normal model example - courtesy of German Rodriguez Chapter 3 HEexample.R: script to run the Bayesian model to analyse the health economic problem described in the chapter (chemotherapy) and the several cost-effectiveness analyses presented throughout the chapter. This example is used throughout chapter 4 as well model.txt: JAGS code for the basic cost-effectiveness analysis modelEVPPI_rho.txt: JAGS code for the analysis of the Expected Value of Partial Perfect Information (EVPPI) for the parameter ρ modelEVPPI_gamma.txt: JAGS code for the analysis of the Expected Value of Partial Perfect Information (EVPPI) for the parameter γ Chapter 4 modelNormal.R: script to run the analysis of the Normal model (pages 129-141); continues the analysis from chapter 2 modelNormalBlocking.txt: JAGS code to run the model using blocking to improve convergence (page 133) modelNormal2.txt: JAGS code to run the compute the predictive distribution (page 135) Chapter 5 Example 1: RCT of acupuncture for chronic headache in primary care acupuncture.R: script to run the cost-effectiveness analysis of acupuncture. Based on this paper dataRCTacupuncture.csv: Dataset used for the acupuncture example - courtesy of David Wonderling, Richard Nixon and Richard Grieve actptRCT.txt: JAGS code to run the normal/normal (on the logit/log scale) model actptRCT_gamma.txt: JAGS code to run the normal/gamma (on the logit/natural scale) model actptRCT_logN.txt: JAGS code to run the normal/log-normal (on the logit/natural scale) model Example 2: Neuraminidase inhibitors to reduce influenza in healthy adults EvSynth.R: script to run the cost-effectiveness analysis based on evidence synthesis for the influenza treatment with neuraminidase EvSynth.txt: JAGS code to run the evidence synthesis model Example 3: Markov model for the treatment of asthma MarkovModel.R: script to run the cost-effectiveness analysis for the treatment of asthma MarkovModel.txt: JAGS code to run the conjugated Markov model

1. ### Set of utility functions to perform Bayesian estimations and health economic analysis
   
2. 
   
3. ############################################################
   
4. ## Produce a traceplot to assess convergence of a MCMC run
   
5. ## node is a string with the name of the node to be plotted
   
6. ## model is the name of the object containing the MCMC simulations
   
7. ## Copyright Gianluca Baio 2012
   
8. mytraceplot <- function(node,model=m,title="",lab=""){
   
9.         xlab <- "Iteration"
   
10.         cmd <- ifelse(class(model)=="rjags",mdl <- model$BUGSoutput,mdl <- model)
    
11.         col <- colors()
    
12.         if (mdl$n.chains==1) {
    
13.                 plot(mdl$sims.array[,1,node],t="l",col=col[which(col=="blue")],xlab=xlab,
    
14.                         ylab=lab,main=title)
    
15.         }
    
16.         if (mdl$n.chains==2) {
    
17.                 cols <- c("blue","red")
    
18.                 plot(mdl$sims.array[,1,node],t="l",col=col[which(col==cols[1])],xlab=xlab,
    
19.                         ylab=lab,main=title,ylim=range(mdl$sims.array[,1:2,node]))
    
20.                 points(mdl$sims.array[,2,node],t="l",col=col[which(col==cols[2])])
    
21.         }
    
22.         if (mdl$n.chains>2) {
    
23.                 cols <- c("blue","red","green","magenta","orange","brown","azure")
    
24.                 plot(mdl$sims.array[,1,node],t="l",col=col[which(col==cols[1])],xlab=xlab,
    
25.                         ylab=lab,main=title,
    
26.                         ylim=range(mdl$sims.array[,1:mdl$n.chains,node]))
    
27.                 for (i in 2:mdl$n.chains) {
    
28.                         points(mdl$sims.array[,i,node],t="l",col=col[which(col==cols[i])])
    
29.                 }
    
30.         }
    
31. }
    
32. 
    
33. 
    
34. 
    
35. ############################################################
    
36. ## Compute the value of parameters (a,b) for a Beta distribution to have mean and sd (m,s)
    
37. ## Copyright Gianluca Baio 2012
    
38. betaPar <- function(m,s){
    
39.         a <- m*( (m*(1-m)/s^2) -1 )
    
40.         b <- (1-m)*( (m*(1-m)/s^2) -1 )
    
41.         list(a=a,b=b)
    
42. }
    
43. 
    
44. 
    
45. 
    
46. ############################################################
    
47. ## Compute the value of parameters (mulog,sigmalog) for a logNormal distribution to have mean and sd (m,s)
    
48. ## Copyright Gianluca Baio 2012
    
49. lognPar <- function(m,s) {
    
50.         s2 <- s^2
    
51.         mulog <- log(m) - .5*log(1+s2/m^2)
    
52.         s2log <- log(1+(s2/m^2))
    
53.         sigmalog <- sqrt(s2log)
    
54.         list(mulog=mulog,sigmalog=sigmalog)
    
55. }
    
56. 
    
57. 
    
58. 
    
59. ############################################################
    
60. ## Compute the parameters of a Beta distribution, given a prior guess for:
    
61. ##  mode = the mode of the distribution
    
62. ##  upp  = an upper bound value for the distribution
    
63. ##  prob = the estimated probability that (theta <= upp)
    
64. ## Based on "Bayesian ideas and data analysis", page 100.
    
65. ## Optimisation method to identify the values of a,b that give required conditions on the Beta distribution
    
66. ## Copyright Gianluca Baio 2012
    
67. betaPar2 <- function(mode,upp,prob){
    
68. N <- 10000
    
69. b <- 1:N
    
70. a <- (1+mode*(b-2))/(1-mode)
    
71. sim <- qbeta(prob,a,b)
    
72. m <- ifelse(prob>=.5,max(which(sim>=upp)),min(which(sim>=upp)))
    
73. M <- ifelse(prob>=.5,min(which(sim<=upp)),max(which(sim<=upp)))
    
74. 
    
75. b <- min(m,M)+(b/N)
    
76. a <- (1+mode*(b-2))/(1-mode)
    
77. sim <- qbeta(prob,a,b)
    
78. m <- ifelse(prob>=.5,max(which(sim>=upp)),min(which(sim>=upp)))
    
79. M <- ifelse(prob>=.5,min(which(sim<=upp)),max(which(sim<=upp)))
    
80. a <- ifelse(m==M,a[m],mean(a[m],a[M]))
    
81. b <- ifelse(m==M,b[m],mean(b[m],b[M]))
    
82. 
    
83. step <- 0.001
    
84. theta <- seq(0,1,step)
    
85. density <- dbeta(theta,a,b)
    
86. 
    
87. norm.dens <- density/sum(density)
    
88. cdf <- cumsum(norm.dens)
    
89. M <- min(which(cdf>=.5))
    
90. m <- max(which(cdf<=.5))
    
91. 
    
92. theta.mode <- theta[which(density==max(density))]
    
93. theta.mean <- a/(a+b)
    
94. theta.median <- mean(theta[m],theta[M])
    
95. theta.sd <- sqrt((a*b)/(((a+b)^2)*(a+b+1)))
    
96. beta.params <- c(a,b,theta.mode,theta.mean,theta.median,theta.sd)
    
97. res1 <- beta.params[1]
    
98. res2 <- beta.params[2]
    
99. theta.mode <- beta.params[3]
    
100. theta.mean <- beta.params[4]
     
101. theta.median <- beta.params[5]
     
102. theta.sd <- beta.params[6]
     
103. list(
     
104. res1=res1,res2=res2,theta.mode=theta.mode,theta.mean=theta.mean,theta.median=theta.median,theta.sd=theta.sd)
     
105. }
     
106. 
     
107. 
     
108. 
     
109. ############################################################
     
110. ## Posterior summary
     
111. ## Creates summary statistics a' la WinBUGS --- currently works only for vectors
     
112. ## Copyright Gianluca Baio 2012
     
113. stats <- function(x){
     
114. c(mean(x),sd(x),quantile(x,.025),median(x),quantile(x,.975))
     
115. }

复制代码

1. ## PROGRAMS THE GIBBS SAMPLING - pag 64-67
   
2. 
   
3. # Define the number of observations & simulations
   
4. n <- 30
   
5. 
   
6. # Load the data
   
7. #y <- rnorm(n,4,10)        # This would generate random values
   
8. y <- c(1.2697,7.7637,2.2532,3.4557,4.1776,6.4320,-3.6623,7.7567,5.9032,7.2671,
   
9.         -2.3447,8.0160,3.5013,2.8495,0.6467,3.2371,5.8573,-3.3749,4.1507,4.3092,
   
10.         11.7327,2.6174,9.4942,-2.7639,-1.5859,3.6986,2.4544,-0.3294,0.2329,5.2846)
    
11. 
    
12. # Defines the hyper-parameters to build the full conditionals
    
13. ybar <- mean(y)
    
14. mu_0 <- 0
    
15. sigma2_0 <- 10000
    
16. alpha_0 <- 0.01
    
17. beta_0 <- 0.01
    
18. 
    
19. # Initialises the parameters
    
20. mu <- tau <- numeric()
    
21. sigma2 <- 1/tau
    
22. 
    
23. mu[1] <- rnorm(1,0,3)
    
24. tau[1] <- runif(1,0,3)
    
25. sigma2[1] <- 1/tau[1]
    
26. 
    
27. # Gibbs sampling (samples from the full conditionals)
    
28. nsim <- 1000
    
29. for (i in 2:nsim) {
    
30.         sigma_n <- sqrt(1/(1/sigma2_0 + n/sigma2[i-1]))
    
31.         mu_n <- (mu_0/sigma2_0 + n*ybar/sigma2[i-1])*sigma_n^2
    
32.         mu[i] <- rnorm(1,mu_n,sigma_n)
    
33. 
    
34.         alpha_n <- alpha_0+n/2
    
35.         beta_n <- beta_0 + sum((y-mu[i])^2)/2
    
36.         tau[i] <- rgamma(1,alpha_n,beta_n)
    
37.         sigma2[i] <- 1/tau[i]
    
38. }
    
39. 
    
40. ## Creates a bivariate contour, using the package mvtnorm
    
41. require(mvtnorm)
    
42. theta <- c(mean(mu),mean(sqrt(sigma2)))
    
43. sigma <- c(var(mu),var(sqrt(sigma2)))
    
44. rho <- cor(mu,sqrt(sigma2))
    
45. ins <- c(-10,10)
    
46. x1 <- seq(theta[1]-abs(ins[1]),theta[1]+abs(ins[1]),length.out=1000)
    
47. x2 <- seq(theta[2]-abs(ins[2]),theta[2]+abs(ins[2]),length.out=1000)
    
48. all <- expand.grid(x1, x2)
    
49. Sigma<-matrix(c(sigma[1], sigma[1]*sigma[2]*rho, sigma[1]*sigma[2]*rho, sigma[2]), nrow=2, ncol=2)
    
50. f.x <- dmvnorm(all, mean = theta, sigma = Sigma)
    
51. f.x2 <- matrix(f.x, nrow=length(x1), ncol=length(x2))
    
52. 
    
53. ##Plots of the results at different simulation lengths
    
54. lw <- c(1,1.6)
    
55. 
    
56. # First 10 iterations
    
57. plot(mu,sqrt(sigma2),col="white",pch=20,cex=.3,xlim=range(mu),ylim=range(sqrt(sigma2)),xlab=expression(mu),ylab=expression(sigma),
    
58.         main="After 10 iterations")
    
59. for (i in 1:(9)) {
    
60.         points(c(mu[i],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i]),t="l",col="grey60")
    
61.         points(c(mu[i+1],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i+1]),t="l",col="grey60")
    
62. }
    
63. text(jitter(mu[1:10]),jitter(sqrt(sigma2[1:10])),1:10,col="grey20")
    
64. contour(x = x1, y = x2, z = f.x2, nlevels=5,add=T,col="black",drawlabels=FALSE,lwd=lw[1])
    
65. 
    
66. 
    
67. # First 30 iterations
    
68. plot(mu,sqrt(sigma2),col="white",pch=20,cex=.3,xlim=range(mu),ylim=range(sqrt(sigma2)),xlab=expression(mu),ylab=expression(sigma),
    
69.         main="After 30 iterations")
    
70. for (i in 1:29) {
    
71.         points(c(mu[i],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i]),t="l",col="grey60")
    
72.         points(c(mu[i+1],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i+1]),t="l",col="grey60")
    
73. }
    
74. text(jitter(mu[1:30]),jitter(sqrt(sigma2[1:30])),1:30,col="grey20")
    
75. contour(x = x1, y = x2, z = f.x2, nlevels=5,add=T,col="black",drawlabels=FALSE,lwd=lw[1])
    
76. 
    
77. 
    
78. # First 100 iterations
    
79. plot(mu,sqrt(sigma2),col="white",pch=20,cex=.3,xlim=range(mu),ylim=range(sqrt(sigma2)),xlab=expression(mu),ylab=expression(sigma),
    
80.         main="After 100 iterations")
    
81. for (i in 1:99) {
    
82.         points(c(mu[i],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i]),t="l",col="grey60")
    
83.         points(c(mu[i+1],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i+1]),t="l",col="grey60")
    
84. }
    
85. text(jitter(mu[1:100]),jitter(sqrt(sigma2[1:100])),1:100,col="grey20")
    
86. contour(x = x1, y = x2, z = f.x2, nlevels=5,add=T,col="black",drawlabels=FALSE,lwd=lw[1])
    
87. 
    
88. 
    
89. # All 1000 iterations
    
90. plot(mu,sqrt(sigma2),col="grey60",pch=20,cex=.3,xlim=range(mu),ylim=range(sqrt(sigma2)),xlab=expression(mu),ylab=expression(sigma),
    
91.         main="After 1000 iterations")
    
92. for (i in 1:(nsim-1)) {
    
93.         points(c(mu[i],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i]),t="l",col="grey60")
    
94.         points(c(mu[i+1],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i+1]),t="l",col="grey60")
    
95. }
    
96. text(jitter(mu),jitter(sqrt(sigma2)),1:nsim,col="grey28")
    
97. contour(x = x1, y = x2, z = f.x2, nlevels=5,add=T,col="black",drawlabels=FALSE,lwd=lw[2])

复制代码

1. ## ANALYSIS OF THE NORMAL MODEL - pag 69-73
   
2. 
   
3. # define the working directory assuming - the user is already in ~/WebMaterial/ch2
   
4. working.dir <- paste(getwd(),"/",sep="")
   
5. 
   
6. # Load data from .dta format
   
7. library(foreign)
   
8. data <- read.dta("phbirths.dta")
   
9. attach(data)
   
10. 
    
11. # Arrange the data and define relevant variables
    
12. y <- grams
    
13. N <- length(y)
    
14. X <- gestate
    
15. k <- 5000
    
16. 
    
17. ## Run JAGS
    
18. library(R2jags)
    
19. 
    
20. ## 1. Uncentred model
    
21. dataJags <- list("N","y","X","k")                                 # define the data
    
22. filein <- paste(working.dir,"modelNormal.txt",sep="")                # define the model file
    
23. params <- c("alpha","beta","sigma")                                # define the parameters to be monitored
    
24. inits <- function(){                                                # create a function to randomly simulate
    
25.         list(alpha=rnorm(1),beta=rnorm(1),lsigma=rnorm(1))        #   initial values for the required nodes
    
26. }
    
27. 
    
28. n.iter <- 10000                                                        # define the number of iterations
    
29. n.burnin <- 9500                                                # define the burn-in
    
30. n.thin <- floor((n.iter-n.burnin)/500)                                # compute the thinning to have 1000 iterations monitored
    
31. 
    
32. # Launch JAGS and print the summary table
    
33. m.1 <- jags(dataJags, inits, params, model.file=filein,
    
34.         n.chains=2, n.iter, n.burnin, n.thin,
    
35.         DIC=TRUE, working.directory=working.dir, progress.bar="text")
    
36. print(m.1,digits=3,intervals=c(0.025, 0.975))
    
37. 
    
38. 
    
39. # 2. Uncentred model with thinning (run for longer chains)
    
40. n.iter <- 50000                                                        # increase number of iterations to improve convergence
    
41. n.burnin <- 9500                                                # same burn-in
    
42. n.thin <- floor((n.iter-n.burnin)/500)                                # compute thinning, this time > 1
    
43. 
    
44. # Launch JAGS and print the summary table
    
45. m.2 <- jags(dataJags, inits, params, model.file=filein,
    
46.         n.chains=2, n.iter, n.burnin, n.thin,
    
47.         DIC=TRUE, working.directory=working.dir, progress.bar="text")
    
48. print(m.2,digits=3,intervals=c(0.025, 0.975))
    
49. 
    
50. 
    
51. # 3. Centred model
    
52. X <- gestate-mean(gestate)                                        # create centred covariate to improve convergence
    
53. dataJags <- list("N","y","X","k")                                 # need to re-define data list
    
54. 
    
55. n.iter <- 10000                                                        # doesn't need so many iterations, now
    
56. n.burnin <- 9500
    
57. n.thin <- floor((n.iter-n.burnin)/500)                                # thinning computed and set to 1 again
    
58. 
    
59. # Launch JAGS and print the summary table
    
60. m.3 <- jags(dataJags, inits, params, model.file=filein,
    
61.         n.chains=2, n.iter, n.burnin, n.thin,
    
62.         DIC=TRUE, working.directory=working.dir, progress.bar="text")
    
63. print(m.3,digits=3,intervals=c(0.025, 0.975))
    
64. 
    
65. 
    
66. ## Make plots
    
67. source("../Utils.R")                                                # load utility functions
    
68. 
    
69. # 1. Uncentred model
    
70. # Traceplot
    
71. par(mfrow=c(2,1))
    
72. mytraceplot("alpha",title="Uncentred model",lab=expression(alpha),model=m.1)
    
73. mytraceplot("beta",lab=expression(beta),model=m.1)
    
74. 
    
75. # Makes ACF plot
    
76. alpha.u <- c(m.1$BUGSoutput$sims.array[,1,1],m.1$BUGSoutput$sims.array[,2,1])
    
77. beta.u <- c(m.1$BUGSoutput$sims.array[,1,2],m.1$BUGSoutput$sims.array[,2,2])
    
78. txt <- substitute("Autocorrelation function for "*alpha* " - Uncentred model")
    
79. acf(alpha.u,main=txt,ylab="")
    
80. 
    
81. # Computes correlation for this model
    
82. cor(alpha.u,beta.u)
    
83. txt2 <- substitute("Scatterplot for "*alpha* " and "*beta* " - Uncentred model")
    
84. plot(m.1$BUGSoutput$sims.list$alpha,m.1$BUGSoutput$sims.list$beta,xlab=expression(alpha),
    
85.         ylab=expression(beta),main=txt2)
    
86. 
    
87. 
    
88. # 2. Uncentred model with thinning
    
89. # Traceplot
    
90. par(mfrow=c(2,1))
    
91. mytraceplot("alpha",title="Uncentred model with thinning",lab=expression(alpha),model=m.2)
    
92. mytraceplot("beta",lab=expression(beta),model=m.2)
    
93. 
    
94. # Makes ACF plot
    
95. alpha.l <- c(m.2$BUGSoutput$sims.array[,1,1],m.2$BUGSoutput$sims.array[,2,1])
    
96. beta.l <- c(m.2$BUGSoutput$sims.array[,1,2],m.2$BUGSoutput$sims.array[,2,2])
    
97. txt <- substitute("Autocorrelation function for "*alpha* " - Uncentred model with thinning")
    
98. acf(alpha.l,main=txt,ylab="")
    
99. 
    
100. # Computes correlation for this model
     
101. cor(alpha.l,beta.l)
     
102. txt4 <- substitute("Scatterplot for "*alpha* " and "*beta* " - Uncentred model with thinning")
     
103. plot(m.2$BUGSoutput$sims.list$alpha,m.2$BUGSoutput$sims.list$beta,xlab=expression(alpha),        
     
104.         ylab=expression(beta),main=txt4)
     
105. 
     
106. 
     
107. # 3. Centred model
     
108. par(mfrow=c(2,1))
     
109. mytraceplot("alpha",model=m.3,title="Centred model",lab=expression(alpha))
     
110. mytraceplot("beta",model=m.3,lab=expression(beta))
     
111. 
     
112. # Makes ACF plot
     
113. alpha.c <- c(m.3$BUGSoutput$sims.array[,1,1],m.3$BUGSoutput$sims.array[,2,1])
     
114. beta.c <- c(m.3$BUGSoutput$sims.array[,1,2],m.3$BUGSoutput$sims.array[,2,2])
     
115. txt <- substitute("Autocorrelation function for "*alpha* " - Centred model")
     
116. acf(alpha.c,main=txt,ylab="")
     
117. 
     
118. # Plots Scatterplot for correlation
     
119. cor(alpha.c,beta.c)
     
120. txt3 <- substitute("Scatterplot for "*alpha* " and "*beta* " - Centred model")
     
121. plot(m.3$BUGSoutput$sims.list$alpha,m.3$BUGSoutput$sims.list$beta,xlab=expression(alpha),
     
122.         ylab=expression(beta),main=txt3)

复制代码

1. ## HEALTH ECONOMIC EVALUATION --- EXAMPLE THROUGHOUT CHAPTER 3 AND CHAPTER 4
   
2. 
   
3. ## Loosely based on Fox-Rushby & Cairns (2005). Economic Evaluation. Open University Press
   
4. ## Parameters of the model
   
5. # pi = probability of side effects (specific to each treatment)
   
6. # rho = reduction in probability of side effects for t=2
   
7. # phi = probability of ambulatory visit, given side effects (constant across treatments)
   
8. # c.amb = cost of ambulatory visit after side effects
   
9. # c.hosp = cost of hospitalisation after side effects
   
10. # c.drug = cost of each treatment (known constant)
    
11. 
    
12. # define the working directory assuming - the user is already in ~/WebMaterial/ch2
    
13. working.dir <- paste(getwd(),"/",sep="")
    
14. source("http://www.statistica.it/gianluca/BMHE/WebMaterial/Utils.R")                # loads utility functions
    
15. 
    
16. ## Generates data for the number of side effects under t=0
    
17. N.studies <- 5
    
18. sample.size <- 32
    
19. prop.pi <- .24
    
20. n <- rpois(N.studies,sample.size)
    
21. se <- rbinom(N.studies,n,prop.pi)
    
22. 
    
23. ## Generates the data for the number of patients with side effects requiring ambulatory care
    
24. prop.gamma <- .619
    
25. amb <- rbinom(N.studies,se,prop.gamma)
    
26. 
    
27. ## Computes the hyperparameters for the informative prior distributions
    
28. m.pi <- .5
    
29. s.pi <- sqrt(.125)
    
30. a.pi <- betaPar(m.pi,s.pi)$a
    
31. b.pi <- betaPar(m.pi,s.pi)$b
    
32. 
    
33. m.gamma <- .5
    
34. s.gamma <- sqrt(.125)
    
35. a.gamma <- betaPar(m.gamma,s.gamma)$a
    
36. b.gamma <- betaPar(m.gamma,s.gamma)$b
    
37. 
    
38. m.rho <- 0.8
    
39. s.rho <- 0.2
    
40. tau.rho <- 1/s.rho^2
    
41. 
    
42. mu.amb <- 120
    
43. sd.amb <- 20
    
44. m.amb <- lognPar(mu.amb,sd.amb)$mulog
    
45. s.amb <- lognPar(mu.amb,sd.amb)$sigmalog
    
46. tau.amb <- 1/s.amb^2
    
47. 
    
48. mu.hosp <- 5500
    
49. sd.hosp <- 980
    
50. m.hosp <- lognPar(mu.hosp,sd.hosp)$mulog
    
51. s.hosp <- lognPar(mu.hosp,sd.hosp)$sigmalog
    
52. tau.hosp <- 1/s.hosp^2
    
53. 
    
54. c.drug <- c(110,520)
    
55. 
    
56. # Number of patients in the population
    
57. N <- 1000
    
58. 
    
59. 
    
60. ## Run JAGS
    
61. library(R2jags)
    
62. dataJags <- list("se","amb","n","N.studies","a.pi","b.pi","a.gamma",                # data list
    
63.                 "b.gamma","m.amb","tau.amb","m.hosp",
    
64.                 "tau.hosp","m.rho","tau.rho","N")
    
65. filein <- "http://www.statistica.it/gianluca/BMHE/WebMaterial/ch3/model.txt"        # model file
    
66. params <- c("pi","gamma","c.amb","c.hosp","rho","SE","A","H")                        # parameters list
    
67. inits <- function(){                                                                # randomly initialise relevant variables
    
68.         SE=rbinom(2,N,.2)
    
69.         list(pi=c(runif(1),NA),gamma=runif(1),c.amb=rlnorm(1),
    
70.         c.hosp=rlnorm(1),rho=runif(1),SE=SE,A=rbinom(2,SE,.2))
    
71. }
    
72. 
    
73. n.iter <- 20000
    
74. n.burnin <- 9500
    
75. n.thin <- floor((n.iter-n.burnin)/250)
    
76. 
    
77. chemo <- jags(dataJags, inits, params, model.file=filein,
    
78.         n.chains=2, n.iter, n.burnin, n.thin,
    
79.         DIC=TRUE, working.directory=working.dir, progress.bar="text")
    
80. ## NB: If initial values for SE are not compatible with the rest of the model, JAGS stops with an error.
    
81. ##     The solution is to re-launch the whole model until it does work.
    
82. 
    
83. print(chemo,digits=3,intervals=c(0.025, 0.975))
    
84. attach.bugs(chemo$BUGSoutput)
    
85. 
    
86. 
    
87. ## Cost-effectiveness analysis
    
88. # Defines the variables of cost and effectiveness
    
89. e <- c <- matrix(NA,chemo$BUGSoutput$n.sims,2)
    
90. e <- N - SE
    
91. for (t in 1:2) {
    
92.         c[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
    
93. }
    
94. 
    
95. library(BCEA)
    
96. treats <- c("Old Chemotherapy","New Chemotherapy")
    
97. m <- bcea(e=e,c=c,ref=2,interventions=treats,Kmax=50000)
    
98. summary(m)
    
99. 
    
100. 
     
101. 
     
102. ## Computes EVPPI
     
103. # Analysis of Expected Value of Partial Information using Strong & Oakley's or Sadatsafavi et al.'s methods
     
104. inp <- CreateInputs(chemo)        # Create the inputs (parameters & matrix of MCMC simulations for all of them)
     
105. x.so <- evppi("rho",inp$mat,m,n.blocks=50)                # Uses S&O method
     
106. plot(x.so)
     
107. x.sal <- evppi("rho",inp$mat,m,n.seps=2)                # Uses Sadatsafavi et al method
     
108. plot(x.sal)
     
109. # Compare the two methods
     
110. plot(x.so$k,x.so$evi,t="l",lwd=2,xlab="Willingness to pay",ylab="",main="Expected value of partial information")
     
111. points(x.so$k,x.so$evppi,t="l",col="blue")
     
112. points(x.sal$k,x.sal$evppi,t="l",col="red")
     
113. legend("topleft",c("EVPI","EVPPI for rho (S&O)","EVPPI for rho (S et al)"),cex=.7,bty="n",lty=1,
     
114.         col=c("black","blue","red"))
     
115. 
     
116. 
     
117. # Analysis of Expected Value of Partial Information using 2stage approach
     
118. # Main parameters:
     
119. # rho = reduction in side effects
     
120. # gamma = chance of ambulatory care
     
121. rho.sim <- rho
     
122. gamma.sim <- gamma
     
123. 
     
124. detach.bugs()
     
125. inits <- NULL
     
126. 
     
127. # EVPPI for rho
     
128. # Replicates the C/E analysis for a fixed value of rho at each iteration
     
129. Ustar.phi <- matrix(NA,chemo$BUGSoutput$n.sims,length(m$k))
     
130. Rhat <- matrix(NA,chemo$BUGSoutput$n.sims,length(chemo$BUGSoutput$summary[,"Rhat"])-1)
     
131. for (i in 1:chemo$BUGSoutput$n.sims) {
     
132.         rho <- rho.sim[i]       
     
133. 
     
134.         dataJags <- list("se","amb","n","N.studies","a.pi","b.pi","a.gamma","b.gamma","m.amb",
     
135.                 "tau.amb","m.hosp","tau.hosp","N","rho")
     
136.         filein <- "modelEVPPI_rho.txt"
     
137.         params <- c("pi","gamma","c.amb","c.hosp","SE","A","H")
     
138.         n.iter <- 10000
     
139.         n.burnin <- 9750
     
140.         n.thin <- floor((n.iter-n.burnin)/250)
     
141.         chemo.evppi <- jags(dataJags, inits, params, model.file=filein,
     
142.                 n.chains=2, n.iter, n.burnin, n.thin,
     
143.                 DIC=TRUE, working.directory=working.dir, progress.bar="text")
     
144.         Rhat[i,] <- chemo.evppi$BUGSoutput$summary[,"Rhat"]
     
145.         attach.bugs(chemo.evppi$BUGSoutput)
     
146. 
     
147.         # Performs the economic analysis given the simulations for all the
     
148.         # other parameters, conditionally on rho
     
149.         e.temp <- c.temp <- matrix(NA,chemo$BUGSoutput$n.sims,2)
     
150.         e.temp <- N - SE
     
151.         for (t in 1:2) {
     
152.                 c.temp[,t] <- c.drug[Rt]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
     
153.         }
     
154.         m.evppi <- bcea(e=e.temp,c=c.temp,ref=2,interventions=treats,Kmax=50000,Ktable=25000)
     
155. 
     
156.         # Computes the maximum expected utility for that iteration for each value of k
     
157.         Ustar.phi[i,] <- apply(m.evppi$Ustar,2,mean)
     
158.         rm(m.evppi); detach.bugs()
     
159.         print(i)
     
160. }
     
161. 
     
162. # Computes the average value of the maximum expected utility for each value of k
     
163. Vstar.phi <- apply(Ustar.phi,2,mean)
     
164. Umax <- apply(apply(m$U,c(2,3),mean),1,max)
     
165. 
     
166. # Computes the EVPPI for each value of k
     
167. EVPPI_rho <- Vstar.phi - Umax
     
168. 
     
169. # Plot the results
     
170. plot(m$k,m$evi,t="l",xlab="Willingness to pay",ylab="")
     
171. points(m$k,EVPPI_rho,t="l",lty=2)
     
172. text <- c("EVPI",expression(paste("EVPPI for ",rho,sep="")))
     
173. legend("topleft",text,lty=1:2,cex=.9,box.lty=0)
     
174. 
     
175. # Checks for convergence in all (inner) simulations
     
176. par(mfrow=c(4,3))
     
177. for (i in 1:12) {
     
178.         plot(Rhat[,i],t="l",main=rownames(chemo.evppi$BUGSoutput$summary)[i],ylim=c(.99,1.15))
     
179.         abline(h=1.1)
     
180. }
     
181. 
     
182. 
     
183. # EVPI for gamma
     
184. # Replicates the C/E analysis for a fixed value of rho at each iteration
     
185. Ustar.phi <- matrix(NA,chemo$BUGSoutput$n.sims,length(m$k))
     
186. Rhat <- matrix(NA,chemo$BUGSoutput$n.sims,length(chemo$BUGSoutput$summary[,"Rhat"])-1)
     
187. for (i in 1:chemo$BUGSoutput$n.sims) {
     
188.         gamma <- gamma.sim[i]       
     
189. 
     
190.         dataJags <- list("se","amb","n","N.studies","a.pi","b.pi","m.rho","tau.rho","m.amb",
     
191.                 "tau.amb","m.hosp","tau.hosp","N","gamma")
     
192.         filein <- "modelEVPPI_gamma.txt"
     
193.         params <- c("pi","rho","c.amb","c.hosp","SE","A","H")
     
194.         n.iter <- 10000
     
195.         n.burnin <- 9750
     
196.         n.thin <- floor((n.iter-n.burnin)/250)
     
197.         chemo.evppi <- jags(dataJags, inits, params, model.file=filein,
     
198.                 n.chains=2, n.iter, n.burnin, n.thin,
     
199.                 DIC=TRUE, working.directory=working.dir, progress.bar="text")
     
200.         Rhat[i,] <- chemo.evppi$BUGSoutput$summary[,"Rhat"]
     
201.         attach.bugs(chemo.evppi$BUGSoutput)
     
202. 
     
203.         # Performs the economic analysis given the simulations for all the
     
204.         # other parameters, conditionally on rho
     
205.         e.temp <- c.temp <- matrix(NA,chemo$BUGSoutput$n.sims,2)
     
206.         e.temp <- N - SE
     
207.         for (t in 1:2) {
     
208.                 c.temp[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
     
209.         }
     
210.         m.evppi <- bcea(e=e.temp,c=c.temp,ref=2,interventions=treats,Kmax=50000,Ktable=25000)
     
211. 
     
212.         # Computes the maximum expected utility for that iteration for each value of k
     
213.         Ustar.phi[i,] <- apply(m.evppi$Ustar,2,mean)
     
214.         rm(m.evppi); detach.bugs()
     
215.         print(i)
     
216. }
     
217. 
     
218. # Computes the average value of the maximum expected utility for each value of k
     
219. Vstar.phi <- apply(Ustar.phi,2,mean)
     
220. Umax <- apply(apply(m$U,c(2,3),mean),1,max)
     
221. 
     
222. # Computes the EVPPI for each value of k
     
223. EVPPI_gamma <- Vstar.phi - Umax
     
224. 
     
225. # Plot the results
     
226. plot(m$k,m$evi,t="l",xlab="Willingness to pay",ylab="")
     
227. points(m$k,EVPPI_gamma,t="l",lty=2)
     
228. text <- c("EVPI",expression(paste("EVPPI for ",gamma,sep="")))
     
229. legend("topleft",text,lty=1:2,cex=.9,box.lty=0)
     
230. 
     
231. 
     
232. # Checks for convergence in all (inner) simulations
     
233. par(mfrow=c(4,3))
     
234. for (i in 1:12) {
     
235.         plot(Rhat[,i],t="l",main=rownames(chemo.evppi$BUGSoutput$summary)[i],ylim=c(.99,1.15))
     
236.         abline(h=1.1)
     
237. }
     
238. 
     
239. 
     
240. 
     
241. ## Model averaging
     
242. # M1 = complete model
     
243. # M2 = rho =: 0.3
     
244. # M3 = rho =: 1
     
245. 
     
246. # analysis for M2
     
247. detach.bugs()
     
248. rho <- 1
     
249. 
     
250. dataJags <- list("se","n","amb","N.studies","a.pi","b.pi","a.gamma","b.gamma","m.amb","tau.amb","m.hosp","tau.hosp","N","rho")
     
251. filein <- "modelEVPPI.txt"
     
252. params <- c("pi","gamma","c.amb","c.hosp","rho","SE","A","H")
     
253. inits <- function(){
     
254.         SE=rbinom(2,N,.2)
     
255.         list(pi=c(runif(1),NA),gamma=runif(1),c.amb=rlnorm(1),c.hosp=rlnorm(1),SE=SE,A=rbinom(2,SE,.2))
     
256. }
     
257. 
     
258. n.iter <- 20000
     
259. n.burnin <- 9500
     
260. n.thin <- floor((n.iter-n.burnin)/250)
     
261. chemo2 <- jags(dataJags, inits, params, model.file=filein,
     
262.         n.chains=2, n.iter, n.burnin, n.thin,
     
263.         DIC=TRUE, working.directory=working.dir, progress.bar="text")
     
264. print(chemo2,digits=3,intervals=c(0.025, 0.975))
     
265. attach.bugs(chemo2$BUGSoutput,overwrite=TRUE)
     
266. 
     
267. ## Cost-effectiveness analysis
     
268. # Defines the variables of cost and effectiveness
     
269. e <- c <- matrix(NA,chemo$BUGSoutput$n.sims,2)
     
270. e <- N - SE
     
271. for (t in 1:2) {
     
272.         c[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
     
273. }
     
274. 
     
275. m2 <- bcea(e=e,c=c,ref=2,interventions=treats,Kmax=50000)
     
276. fileout <- "Model_rho_fixed.ps"
     
277. postscript(fileout,paper="special",width=6,height=6,horizontal=F)
     
278. plot(m2)
     
279. dev.off()
     
280. 
     
281. ## Compares the models based on the DIC
     
282. d <- w <- numeric()
     
283. d[1] <- chemo$BUGSoutput$DIC
     
284. d[2] <- chemo2$BUGSoutput$DIC
     
285. dmin <- min(d)
     
286. 
     
287. # Computes the model weights
     
288. w <- exp(-.5*(d-dmin))/sum(exp(-.5*(d-dmin)))
     
289. 
     
290. # Computes the model average health economic analysis
     
291. SE <- w[1]*chemo$BUGSoutput$sims.list$SE + w[2]*chemo2$BUGSoutput$sims.list$SE
     
292. c.hosp <- w[1]*chemo$BUGSoutput$sims.list$c.hosp + w[2]*chemo2$BUGSoutput$sims.list$c.hosp
     
293. c.amb <- w[1]*chemo$BUGSoutput$sims.list$c.amb + w[2]*chemo2$BUGSoutput$sims.list$c.amb
     
294. A <- w[1]*chemo$BUGSoutput$sims.list$A + w[2]*chemo2$BUGSoutput$sims.list$A
     
295. H <- w[1]*chemo$BUGSoutput$sims.list$H + w[2]*chemo2$BUGSoutput$sims.list$H
     
296. 
     
297. e <- c <- matrix(NA,chemo$BUGSoutput$n.sims,2)
     
298. e <- N - SE
     
299. for (t in 1:2) {
     
300.         c[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
     
301. }
     
302. m.avg <- bcea(e=e,c=c,ref=2,interventions=treats,Kmax=50000)
     
303. plot(m.avg)
     
304. 
     
305. par(mfrow=c(2,2))
     
306. contour(m)
     
307. contour(m2)
     
308. contour(m.avg)

复制代码

1. ## Health economic evaluation --- example for Chapter 3
   
2. ## Loosely based on Fox-Rushby & Cairns (2005). Economic Evaluation. Open University Press
   
3. 
   
4. model {
   
5.         for (s in 1:N.studies) {
   
6.                 se[s] ~ dbin(pi[1],n[s])
   
7.                 amb[s] ~ dbin(gamma,se[s])
   
8.         }
   
9.         pi[1] ~ dbeta(a.pi,b.pi)
   
10.         pi[2] <- pi[1]*rho
    
11.         rho ~ dnorm(m.rho,tau.rho)
    
12.         gamma ~ dbeta(a.gamma,b.gamma)
    
13.         c.amb ~ dlnorm(m.amb,tau.amb)
    
14.         c.hosp ~ dlnorm(m.hosp,tau.hosp)
    
15. 
    
16. # Predictive distributions on the clinical outcomes
    
17.         for (t in 1:2) {
    
18.                 SE[t] ~ dbin(pi[t],N)
    
19.                 A[t] ~ dbin(gamma,SE[t])
    
20.                 H[t] <- SE[t] - A[t]
    
21.         }
    
22. }

复制代码

1. ## Health economic evaluation --- example for Chapter 3 --- Model for the EVPPI
   
2. ## Loosely based on Fox-Rushby & Cairns (2005). Economic Evaluation. Open University Press
   
3. 
   
4. model {
   
5.         for (s in 1:N.studies) {
   
6.                 se[s] ~ dbin(pi[1],n[s])
   
7.                 amb[s] ~ dbin(gamma,se[s])
   
8.         }
   
9.         pi[1] ~ dbeta(a.pi,b.pi)
   
10.         pi[2] <- pi[1]*rho
    
11.         gamma ~ dbeta(a.gamma,b.gamma)
    
12.         c.amb ~ dlnorm(m.amb,tau.amb)
    
13.         c.hosp ~ dlnorm(m.hosp,tau.hosp)
    
14. 
    
15.         for (t in 1:2) {
    
16.                 SE[t] ~ dbin(pi[t],N)
    
17.                 A[t] ~ dbin(gamma,SE[t])
    
18.                 H[t] <- SE[t] - A[t]
    
19.         }
    
20. }

复制代码

1. ## Health economic evaluation --- example for Chapter 3 --- Model for the EVPPI
   
2. ## Loosely based on Fox-Rushby & Cairns (2005). Economic Evaluation. Open University Press
   
3. 
   
4. model {
   
5.         for (s in 1:N.studies) {
   
6.                 se[s] ~ dbin(pi[1],n[s])
   
7.                 amb[s] ~ dbin(gamma,se[s])
   
8.         }
   
9.         pi[1] ~ dbeta(a.pi,b.pi)
   
10.         pi[2] <- pi[1]*rho
    
11.         rho ~ dnorm(m.rho,tau.rho)
    
12.         c.amb ~ dlnorm(m.amb,tau.amb)
    
13.         c.hosp ~ dlnorm(m.hosp,tau.hosp)
    
14. 
    
15.         for (t in 1:2) {
    
16.                 SE[t] ~ dbin(pi[t],N)
    
17.                 A[t] ~ dbin(gamma,SE[t])
    
18.                 H[t] <- SE[t] - A[t]
    
19.         }
    
20. }

复制代码

1. ## ANALYSIS OF THE NORMAL MODEL - pag 129-141
   
2. 
   
3. # define the working directory assuming - the user is already in ~/WebMaterial/ch2
   
4. working.dir <- paste(getwd(),"/",sep="")
   
5. 
   
6. # Load data from .dta format
   
7. library(foreign)
   
8. data <- read.dta("../ch2/phbirths.dta")
   
9. attach(data)
   
10. 
    
11. # Arrange the data and define relevant variables
    
12. y <- grams
    
13. N <- length(y)
    
14. X <- gestate
    
15. k <- 5000
    
16. 
    
17. # Run JAGS
    
18. library(R2jags)
    
19. 
    
20. ## 4. Uncentred model with blocking (for chapter 4)
    
21. ## NB Continues the analysis from chapter 2 (check modelNormal.R in the directory /ch2)
    
22. X <- gestate
    
23. m0 <- c(0,0)
    
24. prec <- 1/100000*matrix(c(1,0,0,1),nrow=2,ncol=2)
    
25. dataJags <- list("m0","prec","N","y","X","k")
    
26. filein <- "modelNormalBlocking.txt"
    
27. params <- c("alpha","beta","sigma")
    
28. inits <- function(){
    
29.         list(lsigma=rnorm(1),coef=rnorm(2,0,1))
    
30. }
    
31. 
    
32. n.iter <- 10000
    
33. n.burnin <- 9500
    
34. n.thin <- floor((n.iter-n.burnin)/500)
    
35. m.4 <- jags(dataJags, inits, params, model.file=filein,
    
36.         n.chains=2, n.iter, n.burnin, n.thin,
    
37.         DIC=TRUE, working.directory=working.dir, progress.bar="text")
    
38. print(m.4,digits=3,intervals=c(0.025, 0.975))
    
39. 
    
40. 
    
41. 
    
42. ## 5. Predictive distribution (for chapter 4)
    
43. # 5a. Uncentered version with thinning
    
44. X.star <- 28
    
45. dataJags2 <- list("N","y","X","k","X.star")
    
46. filein <- paste(working.dir,"modelNormal2.txt",sep="")
    
47. params <- c("alpha","beta","sigma","y.star")
    
48. inits <- function(){
    
49.         list(alpha=rnorm(1),beta=rnorm(1),lsigma=rnorm(1),y.star=runif(1))
    
50. }
    
51. 
    
52. n.iter <- 50000
    
53. n.burnin <- 9500
    
54. n.thin <- floor((n.iter-n.burnin)/500)
    
55. m.5 <- jags(dataJags2, inits, params, model.file=filein,
    
56.         n.chains=2, n.iter, n.burnin, n.thin,
    
57.         DIC=TRUE, working.directory=working.dir, progress.bar="text")
    
58. print(m.5,digits=3,intervals=c(0.025, 0.975))
    
59. 
    
60. # 5b. Centred model with missing data
    
61. y <- c(y,NA)
    
62. N <- length(y)
    
63. X <- c(X,28)
    
64. # still call centred variable X so can still use same model file
    
65. #Xbar <- mean(X)
    
66. #X <- X - Xbar
    
67. dataJags3 <- list("N","y","X","k")
    
68. filein <- "../ch2/modelNormal.txt"
    
69. params3 <- c("alpha","beta","sigma","y[1116]")
    
70. inits <- function(){
    
71. list(alpha=rnorm(1),beta=rnorm(1),lsigma=rnorm(1),
    
72.   y=c(rep(NA,(N-1)),rnorm(1)))
    
73. }
    
74. n.iter <- 50000
    
75. n.burnin <- 9500
    
76. n.thin <- floor((n.iter-n.burnin)/500)
    
77. 
    
78. m.6 <- jags(dataJags3, inits, params3, model.file=filein,
    
79.         n.chains=2,n.iter=n.iter,n.burnin=n.burnin,n.thin=n.thin,DIC=TRUE)
    
80. print(m.6,digits=3,intervals=c(0.025, 0.975))

复制代码