# program to calculate the standard error on the genetic variance in Ve, i.e. varav

varav_exp<-0.05	# genetic variance Ve exponential model
h2m<-0.3		# heritability phenotype
m<-100		# number of families
n<-100		# family size
aw<-0.25	# additive genetic relationship among family members (0.25 for halfsibs; 0.5 for fullsibs; 1.0 for clones)

#derived parameters
t<-aw*h2m
Ve_varW<-(1-h2m)/(1-t)

se_gamma2<-(sqrt(8/m))/n
gamma_squared<-aw*varav*(Ve_varW^2)

se_varav<-se_gamma2*1/aw*(1/Ve_varW^2)

# the optimal family size
optimfam<-2/gamma_squared

# question 1
n=0
results_se=data.frame()

for (i in 1:20) {
n=n+10
t<-aw*h2m
Ve_varW<-(1-h2m)/(1-t)

se_gamma2<-(sqrt(8/m))/n
gamma_squared<-aw*varav*(Ve_varW^2)

se_varav<-se_gamma2*1/aw*(1/Ve_varW^2)
results_se <- rbind(results_se,data.frame(famsize=n,se_varav=se_varav ))
}

directory=getwd()
write("OUTPUT: results_se written to file","results_se.txt",append=FALSE)
write(c(directory,date()),"results_se.txt",append=TRUE)
write.table(results_se,"results_se.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_se.txt written to directory ",directory))


# question 3
varav=0.05
h2m<-0.3		# heritability phenotype
m<-100		# number of families
n<-100		# family size
aw<-0.25	# additive genetic relationship among family members (0.25 for halfsibs; 0.5 for fullsibs; 1.0 for clones)

#derived parameters

h2m=0
results_se=data.frame()

for (i in 1:10) {
h2m=h2m+0.05

t<-aw*h2m
Ve_varW<-(1-h2m)/(1-t)

se_gamma2<-(sqrt(8/m))/n
gamma_squared<-aw*varav*(Ve_varW^2)

se_varav<-se_gamma2*1/aw*(1/Ve_varW^2)
results_se <- rbind(results_se,data.frame(h2m=h2m,se_varav=se_varav ))
}

directory=getwd()
write("OUTPUT: results_se written to file","results_se.txt",append=FALSE)
write(c(directory,date()),"results_se.txt",append=TRUE)
write.table(results_se,"results_se.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_se.txt written to directory ",directory))

# question 4
varav=0.05
h2m<-0.3		# heritability phenotype
m<-100		# number of families
n<-100		# family size
aw<-0.25	# additive genetic relationship among family members (0.25 for halfsibs; 0.5 for fullsibs; 1.0 for clones)

#derived parameters

results_se=data.frame()
m=0
for (i in 1:20) {
m=m+10

t<-aw*h2m
Ve_varW<-(1-h2m)/(1-t)

se_gamma2<-(sqrt(8/m))/n
gamma_squared<-aw*varav*(Ve_varW^2)

se_varav<-se_gamma2*1/aw*(1/Ve_varW^2)
results_se <- rbind(results_se,data.frame(Numfam=m,se_varav=se_varav ))
}

directory=getwd()
write("OUTPUT: results_se written to file","results_se.txt",append=FALSE)
write(c(directory,date()),"results_se.txt",append=TRUE)
write.table(results_se,"results_se.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_se.txt written to directory ",directory))

# question 6
varav=0.05
h2m<-0.3		# heritability phenotype
m<-100		# number of families
n<-100		# family size
aw<-1.0	# additive genetic relationship among family members (0.25 for halfsibs; 0.5 for fullsibs; 1.0 for clones)

#derived parameters

results_se=data.frame()
varav=0
for (i in 1:20) {
varav=varav+0.01

t<-aw*h2m
Ve_varW<-(1-h2m)/(1-t)

se_gamma2<-(sqrt(8/m))/n
gamma_squared<-aw*varav*(Ve_varW^2)

se_varav<-se_gamma2*1/aw*(1/Ve_varW^2)
optimfam<-2/gamma_squared
results_se <- rbind(results_se,data.frame(varav=varav,optimfam=optimfam ))
}

directory=getwd()
write("OUTPUT: results_se written to file","results_se.txt",append=FALSE)
write(c(directory,date()),"results_se.txt",append=TRUE)
write.table(results_se,"results_se.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_se.txt written to directory ",directory))

#question 6 h2m
varav=0.05
h2m<-0.3		# heritability phenotype
m<-100		# number of families
n<-100		# family size
aw<-1.0	# additive genetic relationship among family members (0.25 for halfsibs; 0.5 for fullsibs; 1.0 for clones)

#derived parameters

results_se=data.frame()
h2m=0
for (i in 1:10) {
h2m=h2m+0.05

t<-aw*h2m
Ve_varW<-(1-h2m)/(1-t)

se_gamma2<-(sqrt(8/m))/n
gamma_squared<-aw*varav*(Ve_varW^2)

se_varav<-se_gamma2*1/aw*(1/Ve_varW^2)
optimfam<-2/gamma_squared
results_se <- rbind(results_se,data.frame(h2m=h2m,optimfam=optimfam ))
}

directory=getwd()
write("OUTPUT: results_se written to file","results_se.txt",append=FALSE)
write(c(directory,date()),"results_se.txt",append=TRUE)
write.table(results_se,"results_se.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_se.txt written to directory ",directory))