cat("\n\nTYPE IN AN INTEGER SEED [1-9999999] FOLLOWED BY \"ENTER\"\n")
test<-scan(file="",what=integer(0),n=1)

seed.val<-test
if(seed.val>0)

{
set.seed(seed.val)
}

cat("\n\nTYPE IN THE PREVALENCE OF EXPOSURE [0.001-0.999] FOLLOWED BY \"ENTER\"\n")
test<-scan(file="",what=double(0),n=1)

prevalence.exp<-test

cat("\n\nTYPE IN THE R2 CORRESPONDING TO THE LD [0.001-0.999] FOLLOWED BY \"ENTER\"\n")
test<-scan(file="",what=double(0),n=1)

cat("\n\nPlease wait\n")
R2.target<-test

maxiter<-7
convergence.criterion<-0.0005

#start small and use larger sample sizes near to optimization
sample.size.vector<-c(500000,500000,500000,1000000,2000000)

start.vector<-c(0.00001,0.75,0.9,0.95,0.99,0.999,0.99999)

start.length<-length(start.vector)
midpoint<-round((start.length/2)+1)

reliability<-start.vector

R2<-rep(NA,start.length)
store.reliability<-matrix(NA,nrow=start.length,ncol=(maxiter+1))

store.reliability[,1]<-reliability

convergence<-1
iter<-0

while((abs(convergence)>convergence.criterion)&(iter<maxiter))
{
iter<-iter+1
choose.sample.size<-iter
if(iter>5)choose.sample.size<-5
sample.size<-sample.size.vector[choose.sample.size]
cat("\n\nSTART OF INTERATION",iter,"\n")
cat("Sample size this iteration =",sample.size,"\nSpeed slows with increasing sample size\n\n")


for(j in 1:start.length)
{
print(c(iter,j))
var.b<-reliability[j]
var.w<-(1-var.b)


pheno.orig<-rnorm(sample.size,0,sqrt(var.b))
error.obs1<-rnorm(sample.size,0,sqrt(var.w))
pheno.obs1<-pheno.orig+error.obs1
pheno.true<-pheno.orig

affected.threshold.true<-quantile(pheno.true,prevalence.exp)
affected.threshold.obs1<-quantile(pheno.obs1,prevalence.exp)

affected.true<-(pheno.true<affected.threshold.true)*1
affected.obs1<-(pheno.obs1<affected.threshold.obs1)+1

true.0<-table(affected.true)[1]
true.1<-table(affected.true)[2]

obs1.0<-table(affected.obs1)[1]
obs1.1<-table(affected.obs1)[2]

tab.true.obs1<-table(affected.true,affected.obs1)

sensitivity<-tab.true.obs1[2,2]/true.1
specificity<-tab.true.obs1[1,1]/true.0
R2[j]<-(cor.test(affected.true,affected.obs1)$estimate)^2

if(j==midpoint)
{
pheno.orig.mid<-pheno.orig
error.obs1.mid<-error.obs1
pheno.obs1.mid<-pheno.obs1
error.obs2.mid<-rnorm(sample.size,0,sqrt(var.w))
pheno.obs2.mid<-pheno.orig+error.obs2.mid
reliability.mid<-reliability[j]
sensitivity.mid<-sensitivity
specificity.mid<-specificity
R2.mid<-R2[j]
}

}##j

reliability.save<-reliability

for(t in 1:start.length)
{
rev.t<-(start.length-t)+1
if(R2[t]<R2.target)reliability[1]<-reliability.save[t]
if(R2[rev.t]>R2.target)reliability[start.length]<-reliability.save[rev.t]
}
cat("END OF ITERATION",iter,"\n")


increment<-(reliability[start.length]-reliability[1])/(start.length-1)
convergence<-abs(R2.mid-R2.target)
print("Convergence test")
cat("\n",convergence,"\n\n\n")


for(k in 2:(start.length-1))
{
reliability[k]<-reliability[1]+(k-1)*increment
}

store.reliability[,(iter+1)]<-reliability
}###iter

cat("\n\nCovariance\n")
cov.table<-cov(cbind(pheno.obs1.mid,pheno.obs2.mid))
ICC<-cov.table[1,2]
print(ICC)

cat("\n\nReliability\n")
reliability.mid

cat("\n\nSensitivity\n")
sensitivity.mid

cat("\n\nSpecificity\n")
specificity.mid

cat("\n\nR2\n")
R2.mid

cat("Random seed")
seed.val

cat("\n\n SUMMARY\n #######\n")
cat("\n\n Random Seed\n",seed.val)
cat("\n\n Sample size\n",sample.size)
cat("\n\n Simulated reliability\n",reliability.mid)
cat("\n\n Target R2\n",R2.target)
cat("\n\n Empirical R2\n",R2.mid)
cat("\n\n Empirical reliability\n",ICC)
cat("\n\n Required      Required\n Sensitivity   Specificity\n",sensitivity.mid,"   ",specificity.mid,"\n\n")