################################################################################################################ ### Association between the overall carriage prevalence in adults or 5-17y olds and children under five years ## ################################################################################################################ #### Download JAGS #### library (RCurl) jagsExists = url.exists("http://sourceforge.net/projects/mcmc-jags/files/latest/download") # Regular HTTP if(jagsExists) { txt = getURL("http://sourceforge.net/projects/mcmc-jags/files/latest/download")} #### Use jags package #### library(rjags) #### Set your working directory #### setwd("E:/Pneumo carriage meta-regression") # change as appropriate #### Enter study specific data on number of children included in the sample and number of positive samples #### ncs= 500 # Enter number of children <5y included in the study xcs= 350 # Enter numer of carriers <5y in the study (either overall carriage, or VT, or NVT, as appropriate) #### Data (select as appropriate)#### ## 1. Overall carriage prevalence <5y vs older ## # 1.a. <5y olds vs. adults # studynames<-c("Abdullahi", "Adetifa", "Bello-Gonzales", "Chen", "Dagan", "Granat","Greenberg", "Hammitt", "Henriqus Normark", "Hill", "Hussain", "Inostroza","Kaltoft", "Leino", "Lloyds Evans", "Mueller","Regev-Yochay (2012)", "Regev-Yochay (2004)","Reis") xc<-c(198,375,58,25,59,86,147,377,246,621,87,10,340,15,323,81,189,214,33) # all carriers <5y old in studies with adults nc<-c(349,524,84,94,84,172,216,639,611,666,180,55,584,59,414,128,379,404,50) #all <5y olds in studies with adults xa<-c(16,90,7,0,27,12,33,275,2,821,18,2,23,4,18,28,30,52,19) # all adult carriers na<-c(302,356,64,137,174,154,216,2115,123,1471,237,38,109,123,67,195,376,1300,117) # all adults # 1.b. <5y olds vs. 5-17y olds # #studynames<-c("Hill","Lloyd Evans","Adetifa", "Mueller","Abdullahi","Granat","Parry","Sener", "Hussain", #"Dhakal", "Leino", "Lo", "Cekmez","Reis","Chen", "Inostroza", "Reichler") #xc<-c(621,323,375,81,198,86,192,71,77,18,15,75,0,33,25,10,107) # all carriers <5y old in studies with 5-17y olds #nc<-c(660,414,524,128,349,172,389,248,151,79,59,360,125,50,94,55,166) # all <5y olds in studies with 5-17y olds #xa<-c(621,188,63,57,55,45,212,87,19,26,4,20,25,43,18,5,10) # all carriers 5-17y olds #na<-c(735,342,125,196,213,117,522,413,92,120,31,118,375,95,196,16,53) # all 5-17y olds ## 2. Overall carriage prevalence <1y vs older ## # 2.a. <1y olds vs. adults # #studynames<-c("Abdullahi", "Adetifa", "Darboe","Granat","Hill","Mueller","Nunes","Regev-Yochay (2004)","Regev-Yochay (2012)", "Turner", "van Gils") #xc<-c(58,143,143,49,141,43,83,38,43,188,214) # all carriers <1y old in studies with adults #nc<-c(98,193,196,99,145,62,123,90,90,234,319) #all <1y olds in studies with adults #xa<-c(16,90,26,12,821,28,21,52,30,57,67) # all adult carriers #na<-c(302,356,196,154,1471,195,123,1300,376,231,300) # all adults # 2.b. <1y olds vs. 5-17y olds # #studynames<-c("Abdullahi", "Adetifa", "Granat","Hill","Mueller","Reichler") #xc<-c(58,143,49,141,43,16) # all carriers <1y old in studies with 5-17y olds #nc<-c(98,193,99,145,62,25) # all <1y olds in studies with 5-17y olds #xa<-c(55,63,45,621,57,10) # all carriers 5-17y olds #na<-c(213,125,117,735,196,53) # all 5-17y olds ## 3. VT carriage proportion ## # 3.a. <5y olds vs adults #studynames<-c("Hill", "Adetifa", "Darboe", "Mueller", "Turner", "van Gils", "Reis", "Hammitt", "Regev-Yochay (2012)", "Regev-Yochay (2004)", "Hussain") #xc<-c(108,264,97,45,105,115,12,209,93,89,72) #VT carriage in <5y olds #nc<-c(197,375,143,80,188,213,31,377,189,200,87) #all <5y olds #xa<-c(186,52,6,9,16,27,7,78,9,8,10) #VT carriage in adults #na<-c(674,90,26,28,57,67,19,275,30,29,18) #all adults # 3.b. <5y olds vs. 5-17y olds # #studynames<-c("Hill", "Adetifa", "Mueller", "Reis", "Hussain") #xc<-c(108,264,45,12,72)#VT carriage in <5y olds #nc<-c(197,375,80,31,87) # all <5y olds #xa<-c(226,27,20,12,11)#VT carriage in 5-17y olds #na<-c(617,63,57,43,15)#all 5-17y olds ## 4. VT carriage prevalence ## # 4.a. <5y olds vs adults #studynames<-c("Hill", "Adetifa", "Darboe", "Mueller", "Turner", "van Gils", "Reis", "Hammitt", "Regev-Yochay (2012)", "Regev-Yochay (2004)", "Hussain") #xc<-c(108,264,97,45,105,115,12,209,93,89,72) #VT carriage in <5y olds #nc<-c(219,524,196,128,234,319,50,639,379,404,180) #all <5y olds #xa<-c(186,52,6,9,16,27,7,78,9,8,10) #VT carriage in adults #na<-c(1143,356,196,195,231,300,117,2115,376,1300,237) #all adults # 4.b. <5y olds vs. 5-17y olds # #studynames<-c("Hill", "Adetifa", "Mueller", "Reis", "Hussain") #studynames<-c("Mueller","Hill", "Adetifa", "Mueller", "Hussain", "Reis") #xc<-c(45,108,264,72,12)#VT carriage in <5y olds #nc<-c(128,219,524,180,50) # all <5y olds #xa<-c(20,226,27,11,12)#VT carriage in 5-17y olds #na<-c(196,732,125,71,95)#all 5-17y olds ## 5. NVT carriage prevalence ## # 5.a. <5y olds vs. adults # studynames<-c("Hill", "Adetifa", "Darboe", "Mueller", "Turner", "van Gils", "Reis", "Hammitt", "Regev-Yochay (2012)", "Regev-Yochay (2004)", "Hussain") #xc<-c(89,111,62,40,63,98,21,168,96,111,15) #NVT carriage in <5y olds #nc<-c(219,524,196,128,234,319,50,639,379,404,180) #all <5y olds #xa<-c(488,38,20,9,24,40,12,197,21,22,8) #NVT carriage in adults #na<-c(1143,356,196,231,195,300,117,2115,376,1300,237) #all adults # 5.b. <5y olds vs. 5-17y olds # studynames<-c("Hill", "Adetifa", "Mueller", "Reis", "Hussain") #xc<-c(89,111,40,20,15)#NVT carriage in <5y olds #nc<-c(219,524,128,50,180) # all <5y olds #xa<-c(391,36,37,31,4)#NVT carriage in 5-17y olds #na<-c(732,125,57,95,71)#all 5-17y olds #### Model code ##### jcode <- "model{ for(i in 1:N){ xc[i] ~ dbin(pc[i], nc[i]) #number of carriers in children, follow binom distrib ## Likelihood function derived from the data e.g. 25~bin(theta,100) L=25(power(theta))*(1-25)(power(1-theta)) xa[i] ~ dbin(pa[i], na[i]) # number of carriers in older age groups, follow binom distrib pa[i] ~ dnorm(ma[i], pre) # prevalence in older age groups (pa) distributed around the model prevalence (ma) with a precission (pre) ma[i] <- b0 + b1*pc[i] # the mean of the model. ## Prior for the prevalence in children <5y in each study pc[i] ~ dunif(pcl, pcu) # you can provides the limits along with the data } ## Priors for regression coefficients: b0 ~ dunif(b0l, b0u) b1 ~ dunif(b1l, b1u) pre <- 1/pow(sd,2) # to get the standard deviation of the model sd ~ dunif(sdl, sdu) ## Predictions for the specific survey xcs ~ dbin(pcs, ncs) # likelihood in the survey pcs ~ dunif(pcl, pcu) # prevalence in <5y olds es ~ dnorm(0, pre) # generate random error of model pas <- b0 + b1*pcs + es # prevalence in older age group ## Prevalence estimates for prevalence estimates for each 5% increase in prevalence from 5% to 95% p0<-b0 + b1*0 +es p1 <- b0 + b1*0.1 + es p2<- b0 + b1*0.2 + es p3<- b0 + b1*0.3 + es p4<- b0 + b1*0.4 + es p5<- b0 + b1*0.5 + es p6<- b0 + b1*0.6 + es p7<- b0 + b1*0.7 + es p8<- b0 + b1*0.8+ es p9<- b0 + b1*0.9 + es p05<-b0 + b1*0.05 +es p15<- b0 + b1*0.15 + es p25<- b0 + b1*0.25 + es p35<- b0 + b1*0.35 + es p45<- b0 + b1*0.45 + es p55<- b0 + b1*0.55 + es p65<- b0 + b1*0.65 + es p75<- b0 + b1*0.75 + es p85<- b0 + b1*0.85+ es p95<- b0 + b1*0.95 + es p99<- b0 + b1 + es }" #### Code for running JAGS ### jdat <- list(N=length(nc), xc=xc, nc=nc, xa=xa, na=na,xcs=xcs,ncs=ncs, pcl=0, pcu=1, b0l=-0.5, b0u=1, b1l=-5, b1u=5, sdl=0, sdu=0.1) jmod <- jags.model(textConnection(jcode), data=jdat, n.chains=2, n.adapt=1000) update(jmod,5000) jpos <- coda.samples(jmod, c("b0", "b1", "sd", "pcs", "pas","p0", "p1", "p2", "p3", "p4", "p5", "p6", "p7", "p8", "p9", "p99", "p15", "p25", "p35", "p45", "p55", "p65", "p75", "p85", "p95", "p05"), n.iter=100000, thin=5) #### Retrieve results #### summary(jpos) plot(jpos, ask=TRUE) ## Specific prevalence estimates from the model for the study specific data provided jposmat<-as.matrix(jpos) pcs<-jposmat[,25] # prevalence in <5y olds [from study] pas<-jposmat[,24] # prevalence in older age group (either 5-17y old or adult) [prediction] pcs_quant<-quantile(pcs,probs=c(0.025,0.075,0.125,0.25,0.50,0.75,0.875,0.925,0.975)) # median and 50%, 75%, 85% and 95% prediction intervals pas_quant<-quantile(pas,probs=c(0.025,0.075,0.125,0.25,0.50,0.75,0.875,0.925,0.975)) # median and 50%, 75%, 85% and 95% prediction intervals ## Example of a plot ## (adapt as appropriate) y <- t(apply(t(as.matrix(jpos)[,3:23]), 1, quantile, probs=c(0.025,0.075,0.125,0.25,0.50,0.75,0.875,0.925,0.975))) as.vector(y) x<-matrix(seq(0,1,0.05),21,9) as.vector(x) pc=xc/nc pa=xa/na png(file=paste("regression adults vs <5y olds.png"),width=500, height=400) plot(x, y, pch=".", ylim=c(0,1), col="000000", xlab="Carriage prevalence in children <5y", ylab="Carriage prevalence in adults", axes=F) axis (side=2, at=lx<-seq(0,1,0.1), labels=paste(100*lx,"%",sep="")) axis (side=1, at=ly<-c(0, 0.2, 0.4, 0.6, 0.8, 1), labels=paste(100*ly,"%",sep="")) lines(x[,1],y[,1], col="#961E1A",lty=3) # lower 95% prediction interval limit lines(x[,9],y[,9], col="#961E1A", lty=3) # upper 95% prediction interval limit lines(x[,3],y[,3], col="#08D108", lty=3) # lower 75% prediction interval limit lines(x[,7],y[,7], col="#08D108", lty=3) # upper 75% prediction interval limit lines(x[,8],y[,8], col="#0C22D0", lty=3) # upper 85% prediction interval limit lines(x[,2],y[,2], col="#0C22D0", lty=3) # lower 85% prediction interval limit lines(x[,5],y[,5], col="#423B3B", lty=2) # median points(pc,pa,col="black",pch=1,cex=((na+nc)/max(na+nc))*3+0.5,) dev.off()