### COMBINING Separate estimates of (presumed SAME) parameter
### (as is done in 'fixed-effects' meta-analysis)

# before we combine, lets examine how we would get a single
# estimate, e.g., Redelmeier Study: ADD and car crashes in teenage BOYS

# data from case-control study (BOYS): 
#                     exposed unexposed 
# numerators           767    3421-767   ("case series")
# quasi-denominators   664    3812-664   ("denominator series" or "control series")

# WINBUGS   

model {
#data
cases.exposed       ~ dbin(prop.cases.exposed,3421)  # Binomial, condnl on sum of 2 Poissons
denominaror.exposed ~ dbin(propn.pop.exposed,3812)   # Binomial, by it nature as an independent 'probe'
# structural
prop.cases.exposed <- Rate.Ratio*propn.pop.exposed/(Rate.Ratio*propn.pop.exposed + 1*(1-propn.pop.exposed) )
theta <- log(Rate.Ratio)
# priors
Rate.Ratio ~ dunif(0,10)
propn.pop.exposed  ~ dunif(0,1)	
}

#data 
list(cases.exposed=767, denominaror.exposed=664)

#INITS 
list( Rate.Ratio = 2,  propn.pop.exposed = .05 )


###   R  ####

### COMBINING Separate estimates of (presumed SAME) parameter
### (as is done in 'fixed-effects' meta-analysis)

### Example  ( theta = log[RATE RATIO] )

# Redelmeier Study: ADD and car crashes in teenagers

# BOYS  rate ratio = 1.37, 95% confidence interval 1.22-1.54,

log.rr.hat.boys = (1/100)*round(100*log(1.37),0); log.rr.hat.boys  
se.log.rr.hat.boys = round(log(1.54/1.37)/1.96,2) ; se.log.rr.hat.boys
I.log.rr.hat.boys  = round((1/se.log.rr.hat.boys)^2,0); I.log.rr.hat.boys

# GIRLS rate ratio = 1.31, 95% confidence interval 1.07-1.61

log.rr.hat.girls = (1/100)*round(100*log(1.31),0); log.rr.hat.girls  
se.log.rr.hat.girls = round(log(1.61/1.31)/1.96,2) ; se.log.rr.hat.girls
I.log.rr.hat.girls  = round((1/se.log.rr.hat.girls)^2,0); I.log.rr.hat.girls

# Combine using information-weighting

log.rr.hat=(I.log.rr.hat.boys*log.rr.hat.boys + I.log.rr.hat.girls*log.rr.hat.girls)/
           (I.log.rr.hat.boys                 + I.log.rr.hat.girls)
round(log.rr.hat,3)


# WINBUGS  
 

###### approach A: 

         use 1st point estimate and CI as 'prior' with this amount of Info. 
         add 2nd point estimate [and its Info] as 'data'

model {
theta ~ dnorm (0.31,278)
y ~ dnorm (theta,83)
Rate.Ratio <- exp(theta)
}

#data
list(y=0.27)

#inits
list( theta = -1,0,1)


###### approach B: 

#        use vague prior
#        use 1st point estimate [and its Info] as 'data' 
#        add 2nd point estimate [and its Info] as 'data'

model {
theta ~ dnorm(0,      0.001)
y[1]  ~ dnorm(theta,278    )
y[2]  ~ dnorm(theta,83     )
Rate.Ratio <- exp(theta)
}


#data

list( y=c(0.31, 0.27) )


#inits

list( theta = -10,0,10 )


#### R
#Check Info in Combined estimate (versus 278 + 83)
round(1/fit$summary[1,2]^2,0)



# $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  mac $$$$$$$$$$$

#MCMC via JAGS

MODEL=readLines(,11)   
model {
#data
cases.exposed       ~ dbin(prop.cases.exposed,3421)  # Binomial, condnl on sum of 2 Poissons
denominaror.exposed ~ dbin(propn.pop.exposed,3812)   # Binomial, by it nature as an independent 'probe'
# structural
prop.cases.exposed <- Rate.Ratio*propn.pop.exposed/(Rate.Ratio*propn.pop.exposed + 1*(1-propn.pop.exposed) )
theta <- log(Rate.Ratio)
# priors
Rate.Ratio ~ dunif(0,10)
propn.pop.exposed  ~ dunif(0,1)	
}

writeLines(MODEL,'RateRatioModelBoys.txt')

#data
 
DATA=list(cases.exposed=767, denominaror.exposed=664)

library('rjags')

jags <- jags.model('RateRatioModelBoys.txt',
                   data = DATA,
                   inits=list( Rate.Ratio = 2,  propn.pop.exposed = .05 ),
                   n.chains = 1,
                   n.adapt = 5000)
 
update(jags, 10000)


results=jags.samples(jags, 
  c('Rate.Ratio'), 10000)
str(results)

summary(results$Rate.Ratio[1,,1]) 
hist(results$Rate.Ratio[1,,1],breaks=seq(1.0,1.8,0.02)) 

round( ( sort(results$Rate.Ratio[1,,1]) ) [10000*c(0.025,0.975)] , 2)



#### Combine using MCMC via JAGS  
 

###### approach A: 

         use 1st point estimate and CI as 'prior' with this amount of Info. 
         add 2nd point estimate [and its Info] as 'data'

MODEL=readLines(,5)   
model {
theta ~ dnorm (0.31,278)
y ~ dnorm (theta,83)
Rate.Ratio <- exp(theta)
}

writeLines(MODEL,'RateRatioModelABoysGirls.txt')

#data
DATA=list(y=0.27)

jags <- jags.model('RateRatioModelABoysGirls.txt',
                   data = DATA,
                   inits=list( theta = 2),
                   n.chains = 1,
                   n.adapt = 5000)
 
update(jags, 10000)


results=jags.samples(jags, 
  c('Rate.Ratio'), 10000)
str(results)

summary(results$Rate.Ratio[1,,1]) 

round( ( sort(results$Rate.Ratio[1,,1]) ) [10000*c(0.025,0.500, 0.975)] , 2)


###### approach B: 

#        use vague prior
#        use 1st point estimate [and its Info] as 'data' 
#        add 2nd point estimate [and its Info] as 'data'



MODEL=readLines(,6)   
model {
theta ~ dnorm(0,      0.001)
y[1]  ~ dnorm(theta,278    )
y[2]  ~ dnorm(theta,83     )
Rate.Ratio <- exp(theta)
}

writeLines(MODEL,'RateRatioModelBBoysGirls.txt')


#data

DATA = list( y=c(0.31, 0.27) )

jags <- jags.model('RateRatioModelBBoysGirls.txt',
                   data = DATA,
                   inits=list( theta = 2),
                   n.chains = 1,
                   n.adapt = 5000)
 
update(jags, 10000)


results=jags.samples(jags, 
  c('Rate.Ratio'), 10000)
str(results)

summary(results$Rate.Ratio[1,,1]) 

round( ( sort(results$Rate.Ratio[1,,1]) ) [10000*c(0.025,0.500, 0.975)] , 2)

#### CHECK via logistic regression

#boys

ADD = c(1,0)
case=c(767, 664) ; base=c(2654, 3148)
logistic.fit = glm( cbind(case,base) ~ ADD, family=binomial )

summary(logistic.fit)
round(exp(logistic.fit$coefficients),2)

summary(logistic.fit)$cov.scaled

1/ ( summary(logistic.fit)$cov.scaled)[2,2]  # I = inverse of Var-covar matrix

#Check Info in Combined estimate (versus 278 + 83)