# # HW study in several centers 


kresli <- function(n){
	nopak <- 100;
	tauhat <- numeric(nopak);	
	

	K <- 3;
	tau <- log(2);
#   tau <- 0; 
	
	treatment <- c(rep(1, K*n), rep(0, K*n));
	center <- factor(rep(rep(1:n, rep(K, n)),2));
	
	tau0 <- tau*treatment;
	
	link <- function(x) exp(x)/(1+exp(x));
	
	for(i in 1:nopak){
    psi <- 4*(runif(n) - 0.5);
		psi0 <- rep(rep(psi, rep(K, n)),2);  
# 		browser();
		
		psti <- link(tau0 + psi0)		
		
		X <- (runif(2*K*n) < psti);
# 		X1 <- (runif(K*n) < psti1);
# 		X0 <- (runif(K*n) < psti0);

# 		browser();
		tauhat[i] <- coef(glm(X~treatment+center, family=binomial))[2]
	}
	
# 	browser();
	hist(tauhat, main = paste("I =", n));
	abline(v=tau, col="blue");
	points(tau, 0, col="blue", pch=8)
	points(median(tauhat), 0, col="red", pch=8)
	print(summary(tauhat));
}

pdf(file="Stratified-study-binomial.pdf")
par(mfrow=c(2,2));
kresli(10);
kresli(25);
kresli(50);
kresli(100);
par(mfrow=c(1,1));
dev.off();

# stop();

# # pdf(file="Stratified-study-binomial.pdf")
# par(mfrow=c(2,2));
# kresli(50);
# kresli(100);
# kresli(250);
# kresli(500);
# par(mfrow=c(1,1));
# # dev.off();