rm(list=ls())
source("hap.vc.fun.r")

## (1) get the example datasets

Y.bt  = unlist(read.table("myy.bt", header=T))
Y.qt  = unlist(read.table("myy.qt", header=T))
Geno  = read.table("mygeno", header=T)
X.adj = read.table("myx.adj", header=T)

## (2) Explanation of the arguments needed in the main function "vc.score.fun":

    ##------------------------------------------------------------------------------------------------------
    ## ARGUMENTS of vc.score.fun
    ##------------------------------------------------------------------------------------------------------
    ## y             = vector of trait values
    ## geno          = matrix of alleles in the same format of "haplo.score". That is, each SNP has 2 adjacent
    ##                 columns of alleles, and the order of columns corresponds to the order of SNPs on a chromosome
    ## x.adj         = matrix of non-genetic covariates, EXCLUDING the intercept term
    ## trait.type    = character string defining the type of trait, with values of "gaussian" or "binomial"
    ## approx.method = a number specifying the type of method used to approximate p-value, with values of
    ##                 1=resampling method, 2=Gamma approximation (i.e., the 2-moment approximation), and
    ##                 3=the 3-moment approximation
    ## missing.rm    = T or F. Option "T" performs complete data analysi, i.e., removing individuals if missing
    ##                 in genotypes, trait or covariates. Option "F" imputes the similarity score based on allele
    ##                 frequencies, but still remove individuals with missing traits or missing covariates
    ## n.resamp      = this specifies the resampling sample size for approx.method=1.
    ##------------------------------------------------------------------------------------------------------
    ## OTHER NOTEs:
    ##------------------------------------------------------------------------------------------------------
    ## This program assumes that missing values are coded as "NA"
    ##------------------------------------------------------------------------------------------------------


## (3) Apply the R code on the example datasets:

##----------------------------------------------------------
## Binary trait, imputing missing genotype similarity score
##----------------------------------------------------------
result = vc.score.fun(y=Y.bt, geno=Geno, x.adj=X.adj, trait.type="binomial", approx.method=1, missing.rm=F, n.resamp=1e5)
##      Tvc     pval 
## 21.71929  0.40543 <--- this value should be a little bit off as it is obtained from simulation
result = vc.score.fun(y=Y.bt, geno=Geno, x.adj=X.adj, trait.type="binomial", approx.method=2, missing.rm=F)
##        Tvc       pval 
## 21.7192919  0.3908392 
result = vc.score.fun(y=Y.bt, geno=Geno, x.adj=X.adj, trait.type="binomial", approx.method=3, missing.rm=F)
##        Tvc       pval 
## 21.7192919  0.3722678 


##----------------------------------------------------------
## Binary trait, complete data analysis
##----------------------------------------------------------
result = vc.score.fun(y=Y.bt, geno=Geno, x.adj=X.adj, trait.type="binomial", approx.method=1, missing.rm=T, n.resamp=1e5)
##      Tvc     pval 
## 21.39037  0.38660  <--- this value should be a little bit off as it is obtained from simulation
result = vc.score.fun(y=Y.bt, geno=Geno, x.adj=X.adj, trait.type="binomial", approx.method=2, missing.rm=T)
##        Tvc       pval 
## 21.3903727  0.3940161 
result = vc.score.fun(y=Y.bt, geno=Geno, x.adj=X.adj, trait.type="binomial", approx.method=3, missing.rm=T)
##        Tvc       pval 
## 21.3903727  0.3734369 

##----------------------------------------------------------
## Quantitative trait, imputing missing genotype similarity score
##----------------------------------------------------------
result = vc.score.fun(y=Y.qt, geno=Geno, x.adj=X.adj, trait.type="gaussian", approx.method=1, missing.rm=F, n.resamp=1e5)
##        Tvc       pval 
## 0.09811867 0.68302000  <--- this value should be a little bit off as it is obtained from simulation
result = vc.score.fun(y=Y.qt, geno=Geno, x.adj=X.adj, trait.type="gaussian", approx.method=2, missing.rm=F)
##        Tvc       pval 
## 0.09811867 0.62420429 
result = vc.score.fun(y=Y.qt, geno=Geno, x.adj=X.adj, trait.type="gaussian", approx.method=3, missing.rm=F)
##        Tvc       pval 
## 0.09811867 0.62109950 

##----------------------------------------------------------
## Quantitative trait, complete data analysis
##----------------------------------------------------------
result = vc.score.fun(y=Y.qt, geno=Geno, x.adj=X.adj, trait.type="gaussian", approx.method=1, missing.rm=T, n.resamp=1e5)
##       Tvc      pval 
## 0.1065657 0.6170900  <--- this value should be a little bit off as it is obtained from simulation
result = vc.score.fun(y=Y.qt, geno=Geno, x.adj=X.adj, trait.type="gaussian", approx.method=2, missing.rm=T)
##       Tvc      pval 
## 0.1065657 0.6019660 
result = vc.score.fun(y=Y.qt, geno=Geno, x.adj=X.adj, trait.type="gaussian", approx.method=3, missing.rm=T)
##       Tvc      pval 
## 0.1065657 0.5947237