##                                                    ##
## -------------------------------------------------- ##
## Universität Lübeck                                 ##
## Institut für Medizinische Biometrie und Statistik  ##
## Direktor: Prof. Dr. rer. nat. Andreas Ziegler      ##
## Autor: Oliver Hädike 2006                          ##
## Modifiziert: Friedrich Pahlke 09/2007              ##
## Version: 1.2
## Web: www.imbs-luebeck.de                           ##
## -------------------------------------------------- ##
##                                                    ##
 

theta<-0.0        # recombination fraction between marker and trait locus
alpha<-0.05       # significance level
n<-1082           # sample size (for power calculation)
beta<-0.2         # type II error (for sample size calculation) 1-beta=power
SigmaEq<-1        # residual variance 
rho<-0.0          # residual correlation
p<-0.2            # allele frequency 
H<-0.2424         # broad sense heritability
MOI<-(0)          # mode of inheritence 0- additive  1- dominant (-1)- recessive


##                                                    ##
## -------------------------------------------------- ##
##                                                    ##

SigmaEqundH_in_SGq<-function(SigmaEq,Herit){
 SigmaGq<-Herit*SigmaEq/(1-Herit)
 SigmaGq
}

recA_in_SGq<-function(a,p){
 SGq<-4*p^2*a^2*(1-p^2)
 SGq
}

addA_in_SGq<-function(a,p){
 SGq<-2*p*(1-p)*a^2
 SGq
}

domA_in_SGq<-function(a,p){
 SGq<-4*(1-p)^2*a^2*(1-(1-p)^2)
 SGq
}

recSGq_in_A<-function(SigmaGq,p){
 a<-sqrt(SigmaGq/(4*p^2-4*p^4))
 a
}

addSGq_in_A<-function(SigmaGq,p){
 a<-sqrt(SigmaGq/(2*p*(1-p)))
 a
}

domSGq_in_A<-function(SigmaGq,p){
 a<-sqrt(SigmaGq/(4*(1-p)^2-4*(1-p)^4))
 a
}

DundP_in_SDq<-function(d,p){
 q<-1-p
 d<-(2*p*q*d)^2
 d
}

SigmaE_und_Rho_in_SE<-function(SigmaE,Rho){
 SE<-2*SigmaE*(1-Rho)
 SE
}


Power<-function(alpha,n,Vary0H0,Vary2H0,Vary0H1,Vary2H1,EH1){
 pow<-pt( ( qt(alpha,n-1)*sqrt(Vary0H0+Vary2H0)-0.5*sqrt(n)*EH1)/(sqrt(Vary0H1+Vary2H1)),n-1  )
 pow
}

Fallzahl<-function(alpha,beta,Vary0H0,Vary2H0,Vary0H1,Vary2H1,EH1){
 n<-(  (qnorm(1-beta)*sqrt(Vary0H1+Vary2H1)-qnorm(alpha)*sqrt(Vary0H0+Vary2H0))/(-0.5*EH1) )^2
 n
}

##                                                    ##
## -------------------------------------------------- ##
##                                                    ##


MakeHomKrank<-function(p,a,d,SE){
q<-1-p
mat<-matrix(0,ncol=5,nrow=9)
mat[1,1]<-p^2
mat[1,2]<-p^2+0.5*p*q
mat[1,3]<-p
mat[1,4]<-SE
mat[1,5]<-3*SE^2

mat[2,1]<-p*q
mat[2,2]<-0.5*p*q
mat[2,3]<-0
mat[2,4]<-(a-d)^2+SE
mat[2,5]<-(a-d)^4+6*(a-d)^2*SE+3*SE^2

mat[4,1]<-p*q
mat[4,2]<-0.5*p*q
mat[4,3]<-0
mat[4,4]<-(a-d)^2+SE
mat[4,5]<-(a-d)^4+6*(a-d)^2*SE+3*SE^2

mat[5,1]<-q^2
mat[5,2]<-q^2+0.5*p*q
mat[5,3]<-q
mat[5,4]<-SE
mat[5,5]<-3*SE^2

mat
}

##                                                    ##
## -------------------------------------------------- ##
##                                                    ##

MakeHet<-function(p,a,d,SE){
q<-1-p
mat<-matrix(0,ncol=5,nrow=9)

mat[1,1]<-0
mat[1,2]<-0.25*p^2
mat[1,3]<-0.5*p
mat[1,4]<-SE
mat[1,5]<-3*SE^2

mat[2,1]<-p^2*0.5
mat[2,2]<-0.5*p*q+0.25*p^2
mat[2,3]<-0
mat[2,4]<-(a-d)^2+SE
mat[2,5]<-(a-d)^4+6*(a-d)^2*SE+3*SE^2

mat[3,1]<-0.5*p*q
mat[3,2]<-0
mat[3,3]<-0
mat[3,4]<-4*a^2+SE
mat[3,5]<-16*a^4+24*a^2+3*SE^4

mat[4,1]<-0.5*p^2
mat[4,2]<-0.5*p*q+0.25*p^2
mat[4,3]<-0
mat[4,4]<-(a-d)^2+SE
mat[4,5]<-(a-d)^4+6*(a-d)^2*SE+3*SE^2

mat[5,1]<-q*p
mat[5,2]<-0.25*q^2+0.25*p^2
mat[5,3]<-0.5
mat[5,4]<-SE
mat[5,5]<-3*SE^2

mat[6,1]<-0.5*q^2
mat[6,2]<-0.5*p*q+0.25*q^2
mat[6,3]<-0
mat[6,4]<-(a+d)^2+SE
mat[6,5]<-(a+d)^4+6*(a+d)^2*SE+3*SE^2

mat[7,1]<-0.5*p*q
mat[7,2]<-0
mat[7,3]<-0
mat[7,4]<-4*a^2+SE
mat[7,5]<-16*a^4+24*a^2+3*SE^4

mat[8,1]<-0.5*q^2
mat[8,2]<-0.5*p*q+0.25*q^2
mat[8,3]<-0
mat[8,4]<-(a+d)^2+SE
mat[8,5]<-(a+d)^4+6*(a+d)^2*SE+3*SE^2

mat[9,1]<-0
mat[9,2]<-0.25*q^2
mat[9,3]<-0.5*q
mat[9,4]<-SE
mat[9,5]<-3*SE^2

mat
}

##                                                    ##
## -------------------------------------------------- ##
##                                                    ##

MakeHomGesund<-function(p,a,d,SE){
q<-1-p
mat<-matrix(0,ncol=5,nrow=9)

mat[5,1]<-p^2
mat[5,2]<-p^2+0.5*p*q
mat[5,3]<-p
mat[5,4]<-SE
mat[5,5]<-3*SE^2

mat[6,1]<-p*q
mat[6,2]<-0.5*p*q
mat[6,3]<-0
mat[6,4]<-(a+d)^2+SE
mat[6,5]<-(a+d)^4+6*(a+d)^2*SE+3*SE^2

mat[8,1]<-p*q
mat[8,2]<-0.5*p*q
mat[8,3]<-0
mat[8,4]<-(a+d)^2+SE
mat[8,5]<-(a+d)^4+6*(a+d)^2*SE+3*SE^2

mat[9,1]<-q^2
mat[9,2]<-q^2+0.5*p*q
mat[9,3]<-q
mat[9,4]<-SE
mat[9,5]<-3*SE^2

mat
}

##                                                    ##
## -------------------------------------------------- ##
##                                                    ##

SG<-SigmaEqundH_in_SGq(SigmaEq,H)   
SE<-SigmaE_und_Rho_in_SE(SigmaEq,rho)

if(MOI==0){
a<-addSGq_in_A(SG,p)
d<-0
}
if(MOI==1){
a<-domSGq_in_A(SG,p)
d<-a
}
if(MOI==(-1)){
a<-recSGq_in_A(SG,p)
d<-(-a)
}

sd<-DundP_in_SDq(d,p)
DatenHomozygotGesund<-MakeHomGesund(p,a,d,SE)
DatenHeterozygot<-MakeHet(p,a,d,SE)
DatenHomozygotKrank<-MakeHomKrank(p,a,d,SE)

###
#insert here again for sppsp
###

## wenn homozygot gesund, dann g(Vater) = -a
P_GT_Vater_homo_gesund<-(1-p)^2 

## wenn heterozygot, dann g(Vater) = d
P_GT_Vater_hetero<-2*p*(1-p) 

## wenn homozygot krank, dann g(Vater) = a
P_GT_Vater_homo_krank<-p^2 


##                                                    ##
## -------------------------------------------------- ##
## jetzt die bedingten Erwartungswerte von y          ##
## gegeben IBD für jeden Genotyp des Vaters           ##
## -------------------------------------------------- ##
##                                                    ##

## -------------------------------------------------- ##
## 1. Genotyp des Vaters sei homozygot gesund
E_y_IBD0_22<-0
for( i in 1:9){
	E_y_IBD0_22<-E_y_IBD0_22+DatenHomozygotGesund[i,1]*DatenHomozygotGesund[i,4]
}

E_y_IBD1_22<-0
for( i in 1:9){
	E_y_IBD1_22<-E_y_IBD1_22+DatenHomozygotGesund[i,2]*DatenHomozygotGesund[i,4]
}

E_y_IBD2_22<-0
for( i in 1:9){
	E_y_IBD2_22<-E_y_IBD2_22+DatenHomozygotGesund[i,3]*DatenHomozygotGesund[i,4]
}

## -------------------------------------------------- ##
## 2. Genotyp des Vaters sei heterozygot

E_y_IBD0_12<-0
for( i in 1:9){
	E_y_IBD0_12<-E_y_IBD0_12+DatenHeterozygot[i,1]*DatenHeterozygot[i,4]
}

E_y_IBD1_12<-0
for( i in 1:9){
	E_y_IBD1_12<-E_y_IBD1_12+DatenHeterozygot[i,2]*DatenHeterozygot[i,4]
}

E_y_IBD2_12<-0
for( i in 1:9){
	E_y_IBD2_12<-E_y_IBD2_12+DatenHeterozygot[i,3]*DatenHeterozygot[i,4]
}
 
## -------------------------------------------------- ##
## 3. Genotyp des Vaters sei homozygot krank

E_y_IBD0_11<-0
for( i in 1:9){
	E_y_IBD0_11<-E_y_IBD0_11+DatenHomozygotKrank[i,1]*DatenHomozygotKrank[i,4]
}

E_y_IBD1_11<-0
for( i in 1:9){
	E_y_IBD1_11<-E_y_IBD1_11+DatenHomozygotKrank[i,2]*DatenHomozygotKrank[i,4]
}

E_y_IBD2_11<-0
for( i in 1:9){
	E_y_IBD2_11<-E_y_IBD2_11+DatenHomozygotKrank[i,3]*DatenHomozygotKrank[i,4]
}

###############################################################################
# ACHTUNG NEU !!! VERSUCH DER BERÜCKSICHTIGUNG DER REKOMBINATIONSFREQUENZ !!! #
###############################################################################

psi<-theta^2+(1-theta)^2
RekombMat<-matrix(0,ncol=3,nrow=3)

RekombMat[1,1]<-psi^2
RekombMat[1,2]<-psi*(1-psi)
RekombMat[1,3]<-(1-psi)^2

RekombMat[2,1]<-2*psi*(1-psi)
RekombMat[2,2]<-1-2*psi+2*psi^2
RekombMat[2,3]<-2*psi*(1-psi)

RekombMat[3,1]<-(1-psi)^2
RekombMat[3,2]<-psi*(1-psi)
RekombMat[3,3]<-psi^2

E_y_IBD0_11_t<-RekombMat[1,1]*E_y_IBD0_11+RekombMat[2,1]*E_y_IBD1_11+RekombMat[3,1]*E_y_IBD2_11
E_y_IBD1_11_t<-RekombMat[1,2]*E_y_IBD0_11+RekombMat[2,2]*E_y_IBD1_11+RekombMat[3,2]*E_y_IBD2_11
E_y_IBD2_11_t<-RekombMat[1,3]*E_y_IBD0_11+RekombMat[2,3]*E_y_IBD1_11+RekombMat[3,3]*E_y_IBD2_11

E_y_IBD0_12_t<-RekombMat[1,1]*E_y_IBD0_12+RekombMat[2,1]*E_y_IBD1_12+RekombMat[3,1]*E_y_IBD2_12
E_y_IBD1_12_t<-RekombMat[1,2]*E_y_IBD0_12+RekombMat[2,2]*E_y_IBD1_12+RekombMat[3,2]*E_y_IBD2_12
E_y_IBD2_12_t<-RekombMat[1,3]*E_y_IBD0_12+RekombMat[2,3]*E_y_IBD1_12+RekombMat[3,3]*E_y_IBD2_12

E_y_IBD0_22_t<-RekombMat[1,1]*E_y_IBD0_22+RekombMat[2,1]*E_y_IBD1_22+RekombMat[3,1]*E_y_IBD2_22
E_y_IBD1_22_t<-RekombMat[1,2]*E_y_IBD0_22+RekombMat[2,2]*E_y_IBD1_22+RekombMat[3,2]*E_y_IBD2_22
E_y_IBD2_22_t<-RekombMat[1,3]*E_y_IBD0_22+RekombMat[2,3]*E_y_IBD1_22+RekombMat[3,3]*E_y_IBD2_22


## -------------------------------------------------- ##
# nun die bedingten Erwartungswerte von y gegeben den IBD-Wert gegeben das Quantil
# d.h. hier werden die bedingten Wahrscheinlichkeiten P(Genotyp | x > Schwellenwert) verwendet

E_y_IBD0_Schwellenwert<-E_y_IBD0_22_t*P_GT_Vater_homo_gesund + E_y_IBD0_12_t*P_GT_Vater_hetero + E_y_IBD0_11_t*P_GT_Vater_homo_krank

E_y_IBD1_Schwellenwert<-E_y_IBD1_22_t*P_GT_Vater_homo_gesund + E_y_IBD1_12_t*P_GT_Vater_hetero + E_y_IBD1_11_t*P_GT_Vater_homo_krank

E_y_IBD2_Schwellenwert<-E_y_IBD2_22_t*P_GT_Vater_homo_gesund + E_y_IBD2_12_t*P_GT_Vater_hetero + E_y_IBD2_11_t*P_GT_Vater_homo_krank

##                                                    ##
## -------------------------------------------------- ##
## jetzt die bedingten Wahrscheinlichkeiten P(GK | IBD ) bestimmen miitels :
## P(GK|IBD) = Summe_(GT von 11,12,22) P(GT | x>Schwelle)*P(GK|GT, IBD)

P_GK_IBD0<-c()
P_GK_IBD1<-c()
P_GK_IBD2<-c()

for( i in 1:9){
P_GK_IBD0[i]<-P_GT_Vater_homo_gesund*DatenHomozygotGesund[i,1] + P_GT_Vater_hetero*DatenHeterozygot[i,1] + P_GT_Vater_homo_krank*DatenHomozygotKrank[i,1]
}

for( i in 1:9){
P_GK_IBD1[i]<-P_GT_Vater_homo_gesund*DatenHomozygotGesund[i,2] + P_GT_Vater_hetero*DatenHeterozygot[i,2] + P_GT_Vater_homo_krank*DatenHomozygotKrank[i,2]
}

for( i in 1:9){
P_GK_IBD2[i]<-P_GT_Vater_homo_gesund*DatenHomozygotGesund[i,3] + P_GT_Vater_hetero*DatenHeterozygot[i,3] + P_GT_Vater_homo_krank*DatenHomozygotKrank[i,3]
}

##                                                    ##
## -------------------------------------------------- ##
## jetzt schliesslich noch die bedingten Varianzen von y gegeben IBD= 0,1,2 berechnen

## -------------------------------------------------- ##
## 1. Varianz von y gegeben IBD = 0

Var_y_IBD0<-0
for(i in 1:9){
  Var_y_IBD0<-Var_y_IBD0+P_GK_IBD0[i]*(DatenHeterozygot[i,5]-2*E_y_IBD0_Schwellenwert*DatenHeterozygot[i,4]+E_y_IBD0_Schwellenwert^2)
}

## -------------------------------------------------- ##
## 2. Varianz von y gegeben IBD = 1

Var_y_IBD1<-0
for(i in 1:9){
  Var_y_IBD1<-Var_y_IBD1+P_GK_IBD1[i]*(DatenHeterozygot[i,5]-2*E_y_IBD1_Schwellenwert*DatenHeterozygot[i,4]+E_y_IBD1_Schwellenwert^2)
}

## -------------------------------------------------- ##
## 3. Varianz von y gegeben IBD = 2

Var_y_IBD2<-0
for(i in 1:9){
  Var_y_IBD2<-Var_y_IBD2+P_GK_IBD2[i]*(DatenHeterozygot[i,5]-2*E_y_IBD2_Schwellenwert*DatenHeterozygot[i,4]+E_y_IBD2_Schwellenwert^2)
}

##                                                    ##
## -------------------------------------------------- ##
##                                                    ##
## jetzt die Varianzen unter H_0 = Theta=0.5 berechnen

Var_y_IBD0_H0<-0.25*Var_y_IBD0 + 0.5*Var_y_IBD1 + 0.25*Var_y_IBD2

Var_y_IBD2_H0<-0.25*Var_y_IBD0 + 0.5*Var_y_IBD1 + 0.25*Var_y_IBD2



### Korrigierte Berechnung der Varianzen von y gegeben IBD = k ----------------------------------------

E_y_IBD_Schwellenwert <- c(E_y_IBD0_Schwellenwert, E_y_IBD1_Schwellenwert, E_y_IBD2_Schwellenwert)

Var_y_IBD0_new <- Var_y_IBD0 + sum( RekombMat[1,] * ( E_y_IBD_Schwellenwert - sum(E_y_IBD_Schwellenwert * RekombMat[1,]) )^2 )	
Var_y_IBD1_new <- Var_y_IBD1 + sum( RekombMat[2,] * ( E_y_IBD_Schwellenwert - sum(E_y_IBD_Schwellenwert * RekombMat[2,]) )^2 )
Var_y_IBD2_new <- Var_y_IBD2 + sum( RekombMat[3,] * ( E_y_IBD_Schwellenwert - sum(E_y_IBD_Schwellenwert * RekombMat[3,]) )^2 )

Var_y_IBD0_H0_new <- 0.25 * Var_y_IBD0_new + 0.5 * Var_y_IBD1_new + 0.25 * Var_y_IBD2_new	
Var_y_IBD2_H0_new <- 0.25 * Var_y_IBD0_new + 0.5 * Var_y_IBD1_new + 0.25 * Var_y_IBD2_new

### ENDE Korrigierte Berechnung der Varianzen von y gegeben IBD = k ------------------------------------


##                                                    ##
## -------------------------------------------------- ##
##                                                    ##
## zu guter letzt noch den Erwartungswert von beta unter H1 :
	
EH1<-E_y_IBD2_Schwellenwert-E_y_IBD0_Schwellenwert

##                                                    ##
## -------------------------------------------------- ##
##                                                    ##
## und letztendlich Power und Fallzahl berechnen

power<-Power(alpha, n, Var_y_IBD0_H0_new, Var_y_IBD2_H0_new, Var_y_IBD0_new, Var_y_IBD2_new, EH1)
fallzahl<-Fallzahl(alpha, beta, Var_y_IBD0_H0_new, Var_y_IBD2_H0_new, Var_y_IBD0_new, Var_y_IBD2_new, EH1)

##                                                    ##
## -------------------------------------------------- ##
##                                                    ##


#Result Power
power

#Result Sample Size
fallzahl