## Replication and Bayesian updating
setwd('/Volumes/zzhang4/Private/Shared/Research/Replication')
## generate data with certain features
N<-10000
phi<-1
y<-rnorm(N, 0.1, phi)  ## this indicate a s.e. for bar(y)=.1

##dput(y, file='ydatapower.txt')
##y<-dget('ydata.txt')
## y<-sort(y)

## conduct 100 sets of analysis
## random sample
n<-100
ci.r<-ci.b<-NULL
mu.r<-mu.b<-NULL
se.r<-se.b<-NULL

mu0<-0
phi0<-100

R<-100
for (i in sample(1:R)){
	## select data
	s.y<-y[(1+(i-1)*n):(i*n)]

	## doing replication
	ci.r.1<-qnorm(c(.025, .975), mean(s.y), sqrt(phi/n))
	ci.r<-rbind(ci.r, ci.r.1)
	mu.r<-c(mu.r, mean(s.y))
	se.r<-c(se.r, sqrt(phi/n))

	## doing updates
	mu1<-(mean(s.y)*n/phi+mu0/phi0)/(n/phi+1/phi0)
	phi1<-1/(n/phi+1/phi0)
	
	ci.b.1<-qnorm(c(.025, .975), mu1, sqrt(phi1))
	ci.b<-rbind(ci.b, ci.b.1)
	mu.b<-c(mu.b, mu1)
	se.b<-c(se.b, sqrt(phi1))
	
	mu0<-mu1
	phi0<-phi1
}

## cumulative meta analysis
ci.cm<-ci.r[1,]
mu.cm<-mu.r[1]
se.cm<-se.r[1]

for (i in 2:R){
	mu1<-(mu.cm[i-1]/(se.cm[i-1]^2)+mu.r[i]/(se.r[i]^2))/(1/(se.cm[i-1]^2)+1/(se.r[i]^2))
	phi1<-1/(1/(se.cm[i-1]^2)+1/(se.r[i]^2))
	mu.cm<-c(mu.cm, mu1)
	se.cm<-c(se.cm, sqrt(phi1))
	ci.cm.1<-qnorm(c(.025, .975), mu1, sqrt(phi1))
	ci.cm<-rbind(ci.cm, ci.cm.1)
}

## cumulative meta analysis with publication bias
index<-ci.r[, 1]>0 | ci.r[, 2]<0
id<-(1:100)[index]

ci.r.s<-ci.r[index,]
mu.r.s<-mu.r[index]
se.r.s<-se.r[index]

ci.cm.s<-ci.r.s[1,]
mu.cm.s<-mu.r.s[1]
se.cm.s<-se.r.s[1]

for (i in 2:nrow(ci.r.s)){
	mu1<-(mu.cm.s[i-1]/(se.cm.s[i-1]^2)+mu.r.s[i]/(se.r.s[i]^2))/(1/(se.cm.s[i-1]^2)+1/(se.r.s[i]^2))
	phi1<-1/(1/(se.cm.s[i-1]^2)+1/(se.r.s[i]^2))
	mu.cm.s<-c(mu.cm.s, mu1)
	se.cm.s<-c(se.cm.s, sqrt(phi1))
	ci.cm.s.1<-qnorm(c(.025, .975), mu1, sqrt(phi1))
	ci.cm.s<-rbind(ci.cm.s, ci.cm.s.1)
}


## plot the confidence intervals
par(mfrow=c(2,1))

## for replication
plot(0,0,ylim=c(-.5, .5), xlim=c(0, 100),type='n', xlab='Replication', ylab='CI')
abline(h=0)

for (i in 1:R){
	if (min(ci.r[i, ])>0 | max(ci.r[i, ])<0){
		lines(c(i,i),ci.r[i,], lty=3, col='red')
		points(i, mu.r[i], col='red', pch=46, cex=4)
	}else{
		lines(c(i,i),ci.r[i,])
		points(i, mu.r[i], pch=46, cex=4)
	}
}


## for Bayesian analysis

plot(0,0,ylim=c(-.5, .5), xlim=c(0, 100),type='n', xlab='Bayesian', ylab='CI')
abline(h=0)

for (i in 1:R){
	if (min(ci.b[i, ])>0 | max(ci.b[i, ])<0){
		lines(c(i,i),ci.b[i,], lty=3, col='red')
		points(i,mu.b[i], col='red', pch=46, cex=4)
	}else{
		lines(c(i,i), ci.b[i,])
		points(i, mu.b[i], pch=46, cex=4)
	}
}

##plot(1:100, se.b)


## plot the confidence intervals
par(mfrow=c(2,1))

## for replication
plot(0,0,ylim=c(-.5, .5), xlim=c(0, 100),type='n', xlab='Replication', ylab='CI')
abline(h=0)

for (i in 1:R){
	if (min(ci.r[i, ])>0 | max(ci.r[i, ])<0){
		lines(ci.r[i,], c(i,i), lty=3, col='red')
		points(mu.r[i], i, col='red', pch=46, cex=4)
	}else{
		lines(ci.r[i,], c(i,i))
		points(mu.r[i], i, pch=46, cex=4)
	}
}


## for Bayesian analysis

plot(0,0,ylim=c(-.5, .5), xlim=c(0, 100),type='n', xlab='Bayesian', ylab='CI')
abline(h=0)

for (i in 1:R){
	if (min(ci.b[i, ])>0 | max(ci.b[i, ])<0){
		lines(ci.b[i,], c(i,i), lty=3, col='red')
		points(mu.b[i], i, col='red', pch=46, cex=4)
	}else{
		lines(ci.b[i,], c(i,i))
		points(mu.b[i], i, pch=46, cex=4)
	}
}

##plot(1:100, se.b)

## for power
## plot the confidence intervals
par(mfrow=c(2,1))

## for replication
plot(0,0,ylim=c(-.4, .7), xlim=c(0, 100),type='n', xlab='Replication', ylab='CI')
abline(h=0)

for (i in 1:R){
	if (min(ci.r[i, ])>0 | max(ci.r[i, ])<0){
		lines(c(i,i),ci.r[i,], lty=3, col='red')
		points(i, mu.r[i], col='red', pch=46, cex=4)
	}else{
		lines(c(i,i),ci.r[i,])
		points(i, mu.r[i], pch=46, cex=4)
	}
}


## for Bayesian analysis

plot(0,0,ylim=c(-.4, .7), xlim=c(0, 100),type='n', xlab='Bayesian', ylab='CI')
abline(h=0)

for (i in 1:R){
	if (min(ci.b[i, ])>0 | max(ci.b[i, ])<0){
		lines(c(i,i),ci.b[i,], lty=3, col='red')
		points(i,mu.b[i], col='red', pch=46, cex=4)
	}else{
		lines(c(i,i), ci.b[i,])
		points(i, mu.b[i], pch=46, cex=4)
	}
}

#################################
## plot the confidence intervals
## effect (.1)
#################################
## plot the confidence intervals
par(mfrow=c(1,4))

## for replication
plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Replication', ylab='CI')
abline(v=c(.1))
abline(v=c(0), lty=3)

for (i in 1:R){
	if (min(ci.r[i, ])>0 | max(ci.r[i, ])<0){
		lines(ci.r[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.r[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.r[i,], c(i,i))
		points(mu.r[i], i, pch=46, cex=4)
	}
}


## for Bayesian analysis

plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Bayesian', ylab='CI')
abline(v=c(.1))
abline(v=c(0), lty=3)

for (i in 1:R){
	if (min(ci.b[i, ])>0 | max(ci.b[i, ])<0){
		lines(ci.b[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.b[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.b[i,], c(i,i))
		points(mu.b[i], i, pch=46, cex=4)
	}
}

## for cumulative Meta-Analysis
plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Meta Analysis', ylab='CI')
abline(v=c(.1))
abline(v=c(0), lty=3)

for (i in 1:R){
	if (min(ci.cm[i, ])>0 | max(ci.cm[i, ])<0){
		lines(ci.cm[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.cm[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.cm[i,], c(i,i))
		points(mu.cm[i], i, pch=46, cex=4)
	}
}

## for cumulative Meta-Analysis with strong publication bias
plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Meta Analysis (bias)', ylab='CI')
abline(v=c(.1))
abline(v=c(0), lty=3)

for (i in 1:nrow(ci.cm.s)){
	if (min(ci.cm.s[i, ])>0 | max(ci.cm.s[i, ])<0){
		lines(ci.cm.s[i,], c(id[i],id[i]), lty=1, lwd=2, col='red')
		points(mu.cm.s[i], id[i], col='red', pch=46, cex=6)
	}else{
		lines(ci.cm.s[i,], c(id[i],id[i]))
		points(mu.cm.s[i], id[i], pch=46, cex=4)
	}
}

#################################
## plot the confidence intervals
## no effect
#################################
par(mfrow=c(1,4))

## for replication
plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Replication', ylab='CI')
abline(v=c(0))

for (i in 1:R){
	if (min(ci.r[i, ])>0 | max(ci.r[i, ])<0){
		lines(ci.r[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.r[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.r[i,], c(i,i))
		points(mu.r[i], i, pch=46, cex=4)
	}
}


## for Bayesian analysis

plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Bayesian', ylab='CI')
abline(v=c(0))

for (i in 1:R){
	if (min(ci.b[i, ])>0 | max(ci.b[i, ])<0){
		lines(ci.b[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.b[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.b[i,], c(i,i))
		points(mu.b[i], i, pch=46, cex=4)
	}
}

## for cumulative Meta-Analysis
plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Meta Analysis', ylab='CI')
abline(v=c(0))

for (i in 1:R){
	if (min(ci.cm[i, ])>0 | max(ci.cm[i, ])<0){
		lines(ci.cm[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.cm[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.cm[i,], c(i,i))
		points(mu.cm[i], i, pch=46, cex=4)
	}
}

## for cumulative Meta-Analysis with strong publication bias
plot(0,0,xlim=c(-.4, .7), ylim=c(0, 100),type='n', xlab='Meta Analysis (bias)', ylab='CI')
abline(v=c(0))

for (i in 1:nrow(ci.cm.s)){
	if (min(ci.cm.s[i, ])>0 | max(ci.cm.s[i, ])<0){
		lines(ci.cm.s[i,], c(id[i],id[i]), lty=1, lwd=2, col='red')
		points(mu.cm.s[i], id[i], col='red', pch=46, cex=6)
	}else{
		lines(ci.cm.s[i,], c(id[i],id[i]))
		points(mu.cm.s[i], id[i], pch=46, cex=4)
	}
}

###################
## self-correction?
###################

par(mfrow=c(1,4))

## for replication
plot(0,0,xlim=c(-3, 3), ylim=c(0, 100),type='n', xlab='Replication', ylab='CI')
abline(v=c(0))

for (i in 1:R){
	if (min(ci.r[i, ])>0 | max(ci.r[i, ])<0){
		lines(ci.r[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.r[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.r[i,], c(i,i))
		points(mu.r[i], i, pch=46, cex=4)
	}
}


## for Bayesian analysis

plot(0,0,xlim=c(-3, 3), ylim=c(0, 100),type='n', xlab='Bayesian', ylab='CI')
abline(v=c(0))

for (i in 1:R){
	if (min(ci.b[i, ])>0 | max(ci.b[i, ])<0){
		lines(ci.b[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.b[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.b[i,], c(i,i))
		points(mu.b[i], i, pch=46, cex=4)
	}
}

## for cumulative Meta-Analysis
plot(0,0,xlim=c(-3, 3), ylim=c(0, 100),type='n', xlab='Meta Analysis', ylab='CI')
abline(v=c(0))

for (i in 1:R){
	if (min(ci.cm[i, ])>0 | max(ci.cm[i, ])<0){
		lines(ci.cm[i,], c(i,i), lty=1, lwd=2, col='red')
		points(mu.cm[i], i, col='red', pch=46, cex=6)
	}else{
		lines(ci.cm[i,], c(i,i))
		points(mu.cm[i], i, pch=46, cex=4)
	}
}

## for cumulative Meta-Analysis with strong publication bias
plot(0,0,xlim=c(-3, 3), ylim=c(0, 100),type='n', xlab='Meta Analysis (bias)', ylab='CI')
abline(v=c(0))

for (i in 1:nrow(ci.cm.s)){
	if (min(ci.cm.s[i, ])>0 | max(ci.cm.s[i, ])<0){
		lines(ci.cm.s[i,], c(id[i],id[i]), lty=1, lwd=2, col='red')
		points(mu.cm.s[i], id[i], col='red', pch=46, cex=6)
	}else{
		lines(ci.cm.s[i,], c(id[i],id[i]))
		points(mu.cm.s[i], id[i], pch=46, cex=4)
	}
}

### Unweighted results
## conduct 100 sets of analysis
## random sample
n<-100
ci.r<-ci.b<-NULL
mu.r<-mu.b<-NULL
se.r<-se.b<-NULL

mu0<-0
phi0<-100

R<-100
for (i in sample(1:R)){
	## select data
	s.y<-y[(1+(i-1)*n):(i*n)]

	## doing replication
	ci.r.1<-qnorm(c(.025, .975), mean(s.y), sqrt(phi/n))
	ci.r<-rbind(ci.r, ci.r.1)
	mu.r<-c(mu.r, mean(s.y))
	se.r<-c(se.r, sqrt(phi/n))

	## doing updates
	mu1<-(mean(s.y)+mu0)/(2)
	phi1<-1/(n/phi+1/phi0)
	
	ci.b.1<-qnorm(c(.025, .975), mu1, sqrt(phi1))
	ci.b<-rbind(ci.b, ci.b.1)
	mu.b<-c(mu.b, mu1)
	se.b<-c(se.b, sqrt(phi1))
	
	mu0<-mu1
	phi0<-phi1
}

## cumulative meta analysis
ci.cm<-ci.r[1,]
mu.cm<-mu.r[1]
se.cm<-se.r[1]

for (i in 2:R){
	mu1<-(mu.cm[i-1]+mu.r[i])/2
	phi1<-1/(1/(se.cm[i-1]^2)+1/(se.r[i]^2))
	mu.cm<-c(mu.cm, mu1)
	se.cm<-c(se.cm, sqrt(phi1))
	ci.cm.1<-qnorm(c(.025, .975), mu1, sqrt(phi1))
	ci.cm<-rbind(ci.cm, ci.cm.1)
}

## cumulative meta analysis with publication bias
index<-ci.r[, 1]>0 | ci.r[, 2]<0
id<-(1:100)[index]

ci.r.s<-ci.r[index,]
mu.r.s<-mu.r[index]
se.r.s<-se.r[index]

ci.cm.s<-ci.r.s[1,]
mu.cm.s<-mu.r.s[1]
se.cm.s<-se.r.s[1]

for (i in 2:nrow(ci.r.s)){
	mu1<-(mu.cm.s[i-1]+mu.r.s[i])/2
	phi1<-1/(1/(se.cm.s[i-1]^2)+1/(se.r.s[i]^2))
	mu.cm.s<-c(mu.cm.s, mu1)
	se.cm.s<-c(se.cm.s, sqrt(phi1))
	ci.cm.s.1<-qnorm(c(.025, .975), mu1, sqrt(phi1))
	ci.cm.s<-rbind(ci.cm.s, ci.cm.s.1)
}