# Examples:

######################################################################
########################  Example 1  ##################################
######################################################################
### EXAMPLE1. Generate data according to a specific model:

## Epistasis model, no main effects (beta1,beta2=0), interaciton effect gamma=0.4:

> f.epistasis=gen.f(alpha=log(0.1), beta1=0, beta2=0, gamma=0.4, x1=0,x2=0)

> f.epistasis # print f

       [,1]  [,2]  [,3]
[1,] 0.091 0.091 0.091
[2,] 0.091 0.091 0.091
[3,] 0.091 0.091 0.130


## n1=1000 cases; n0=1000 controls, penetrance is specified as f.epistasis,
## minor allele freq (MAF) at Locus is p1=0.4, MAF at Locus 2 is p2=0.3:

> data=gen.data(n1=1000, n0=1000,f=f.epistasis, p1=0.4, p2=0.3)

> data # print data, data is in a 'list' format.

$case
     [,1] [,2] [,3]
[1,]  179  132   35
[2,]  240  207   45
[3,]   81   61   20

$control
     [,1] [,2] [,3]
[1,]  189  161   38
[2,]  238  185   40
[3,]   71   67   11

############################################
## partitioning chisq:
# results are 
# total: total chisq; 
# main: main effect test; 
# interaction: interaction test; 
# par1-par4: four 1-df components for the interaction test (sum of these four
# chisq =interaction.test). These chisq are for epistasis model.
# L1 and L2: marginal tests (Pearson chisq for 2x3 table) at the two loci.

> result = pc(data)

>result

      total        main interaction        par1        par2        par3 
 8.35456166  4.59727084  3.75729081  1.41294189  0.00535921  0.72918531 
       par4          L1          L2 
 1.63168941  3.82730639  0.85209173 


# don't show par1-par4 (they are for epistasis model only):

> result[c(1:3, 8:9)]

      total        main interaction          L1          L2 
  8.3545617   4.5972708   3.7572908   3.8273064   0.8520917 

# remark: if you want par1-par4 for threshold interaction model, then need to switch alleles A and a, B and b. 


######################################################################
########################  Example 2  ##################################
######################################################################
## use real data instead of generate random numbers.
## data can be in the format as shown above or can be a 2x9 table

#input data:

> data1=matrix(c(179, 240,81,132,207,61,35,45,20,189,238,71,161,185,67,38,40,11), nrow=2,ncol=9, byrow=T)

> data1

        aa bb Aa bb AA bb aa Bb Aa Bb AA Bb aa BB Aa BB AA BB
case      179   240    81   132   207    61    35    45    20
control   189   238    71   161   185    67    38    40    11

# or in a list format:

> data2=list(case=matrix(c(179, 240,81,132,207,61,35,45,20),3,3 ), control=matrix(c(189,238,71,161,185,67,38,40,11),3,3   ))

> data2
$case
     [,1] [,2] [,3]
[1,]  179  132   35
[2,]  240  207   45
[3,]   81   61   20

$control
     [,1] [,2] [,3]
[1,]  189  161   38
[2,]  238  185   40
[3,]   71   67   11



#partitioning:

> pc(data1)
> pc(data2)

      total        main interaction        par1        par2        par3 
 8.35456166  4.59727084  3.75729081  1.41294189  0.00535921  0.72918531 
       par4          L1          L2 
 1.63168941  3.82730639  0.85209173