vector<double> &,

vector<double> &,

vector<vector<double> > &,

vector<Individual*>&,

{

if (!par::SNP_major)

Ind2SNP();

// Do not allow monomorphic alleles

if (par::min_af==0)

error("Cannot specify --maf 0 when using --T2; set --maf > 0");

// Do not allow completely missing SNPs

if (par::MAX_GENO_MISSING==1)

error("Cannot specify --geno 1 when using --T2; set --geno < 1");

// Are we using sets? If so, construct these now

if (!par::set_test)

error("You need to specify sets (--set option) with --T2");

// Do not allow quantitative traits

if (!par::bt)

error("Cannot specify --T2 with quantitative traits");

// Prune SET (0-sized sets, MAF==0 SNPs, etc)

pS->pruneSets(*this);

int ns = snpset.size();

///////////////////////////////////////////

// Count how many cases, how many controls

int caseN = 0;

int controlN = 0;

for (int i=0; i < n; i )

if (!sample[i]->missing)

if (sample[i]->aff)

caseN ;

else

controlN ;

if ( caseN == 0 ||

controlN == 0 )

error("No cases / no controls for T(2) test");

// Multi-collinearity SNP pruning

setFlags(true);

pS->pruneMC(*this,true,par::vif_threshold);

// Empirical p-values (1 per set)

perm.setTests(ns);

////////////////////////////////

// Set up permutation structure

// (we need to perform this step

// whether or not we also

// subsequently permute)

perm.setPermClusters(*this);

perm.originalOrder();

vector<double> original = calcHotel(true,

perm,

*pS,

caseN, controlN);

////////////////////////////

// If no permutation, then

// leave now

if (!par::permute) return;

//////////////////////

// Begin permutations

bool finished = false;

while(!finished)

{

// Store permuted results

vector<double> pr(ns);

if (par::perm_genedrop)

perm.geneDrop();

else

perm.permuteInCluster();

pr = calcHotel(false,

perm,

*pS,

caseN, controlN);

////////////////////////////////

// Standard permutation counting

finished = perm.update(pr,original);

} // next permutation

if (!par::silent)

cout << "\n\n";

////////////////////

// Display results

ofstream ASC;

string f;

if (par::adaptive_perm) f = par::output_file_name ".T2.perm";

else f = par::output_file_name ".T2.mperm";

ASC.open(f.c_str(),ios::out);

ASC.precision(4);

printLOG("Writing permutation T2 test results to [ " f " ] \n");

ASC << setw(12) << "SET" << " "

<< setw(4)<< "SIZE" << " "

<< setw(12) << "EMP1" << " ";

if (par::adaptive_perm)

ASC << setw(12)<< "NP" << " ";

else

ASC << setw(12)<< "EMP2" << " ";

ASC << "\n";

for (int s=0; s<snpset.size(); s )

{

ASC << setw(12) << setname[s] << " "

<< setw(4) << snpset[s].size() << " ";

ASC << setw(12) << perm.pvalue(s) << " ";

if (par::adaptive_perm)

ASC << setw(12) << perm.reps_done(s) << " ";

else

ASC << setw(12) << perm.max_pvalue(s) << " ";

ASC << "\n";

}

ASC.close();

Perm & perm,

Set & S,

int ncase,

int ncontrol)

{

ofstream ASC;

if (disp)

{

string f = par::output_file_name ".T2";

ASC.open(f.c_str(),ios::out);

ASC.precision(4);

printLOG("Writing T2 test results to [ " f " ] \n");

ASC << setw(12) << "SET" << " "

<< setw(4)<< "SIZE" << " "

<< setw(12)<< "T2" << " "

<< setw(12) << "DF1" << " "

<< setw(12)<< "DF2" << " "

<< setw(12)<< "P_HOTEL" << "\n";

}

// Number of SETs

int ns = pS->snpset.size();

vector<double> T2(ns,0);

// Consider each SET

for (int s=0; s<ns; s )

{

// Adaptive permutation: skip possibly

if (par::adaptive_perm && !perm.snp_test[s]) continue;

// Consider each SNP, coded 1,0 and -1

// Use mean-substitution for missing alleles

int nss0 = snpset[s].size();

//////////////////////////////////////////////

// Create a vector of pointers for SNPs in set

int nss = 0;

vector<CSNP*> pSNP(0);

for (int j=0; j<nss0; j )

{

// include this SNP? (MC considerations)

if ( pS->cur[s][j] )

{

// Add to set list

pSNP.push_back( SNP[snpset[s][j]] );

// Increase the actual number of snps in set

nss ;

}

}

vector<double> mean2(nss,0); // Case mean

vector<double> mean1(nss,0); // Control mean

vector<vector<double> > pooled; // Covariance matrix

///////////////////////////////////////

// Calculate mean and variance (pooled)

// after imputing missing SNPs

calcHotelSetMeanVariance(pSNP,mean1,mean2,pooled,sample,ncase,ncontrol);

///////////////////////////////

// 2. Calculate test statistic

for (int j1=0; j1<nss; j1 )

for (int j2=j1; j2<nss; j2 )

{

pooled[j1][j2] /= (double)(ncase ncontrol-2);

pooled[j1][j2] *= 1/(double)ncase 1/(double)ncontrol;

if (j1!=j2)

pooled[j2][j1] = pooled[j1][j2];

}

// Get inverse of this matix

bool flag = true;

pooled = svd_inverse(pooled,flag);

// Calculate T2 statistic

vector<double> tmp(nss,0);

for (int j1=0; j1<nss; j1 )

for (int j2=0; j2<nss; j2 )

tmp[j1] = ( mean1[j2] - mean2[j2] ) * pooled[j1][j2];

double stat_t2 = 0;

for (int j1=0; j1<nss; j1 )

stat_t2 = tmp[j1] * ( mean1[j1] - mean2[j1] );

// T(2) is distributed as (n-1)p / (n-p) F(p,n-p)

// where n = number of individuals; p = number of variables

// For two-sample T(2), replace n with n1 n2-1

// Make statistic; test against F(nss, ncase ncontrol-1 - nss)

stat_t2 /=

(

(double)((ncase ncontrol-2)*nss)

/ (double)(ncase ncontrol-nss-1)

);

// Asymptotic p-value, use 1-p as test statistic

// for permutation

double T2p = pF(stat_t2,nss, ncase ncontrol-nss-1);

T2[s] = 1 - T2p;

// cout << "\nTest STAT = " << T2[s] << "\n";

if (disp)

{

ASC << setw(12) << setname[s] << " "

<< setw(4) << snpset[s].size() << " "

<< setw(12) << stat_t2<< " "

<< setw(12) << nss << " "

<< setw(12) << ncase ncontrol-nss-1 << " "

<< setw(12) << T2p << "\n";

}

}

return T2;

vector<double> & mean1,

vector<double> & mean2,

vector<vector<double> > & pooled,

vector<Individual*> & sample,

int ncase,

int ncontrol)

{

int nss = mean1.size();

vector<double> mean(nss,0);

vector<int> cnt1(nss,0);

vector<int> cnt2(nss,0);

////////////////////////////

// Iterate over SNPs in SET

vector<CSNP*>::iterator ps = pSNP.begin();

int j=0;

while ( ps != pSNP.end() )

{

///////////////////////////

// Iterate over individuals

vector<Individual*>::iterator gperson = sample.begin();

vector<bool>::iterator i1 = (*ps)->one.begin();

vector<bool>::iterator i2 = (*ps)->two.begin();

int i=0;

while ( gperson != sample.end() )

{

// Permuted self

Individual * pperson = (*gperson)->pperson;

// Affected individuals

if ( ! pperson->missing )

{

if (pperson->aff)

{

if ( *i1 )

{

if ( *i2 ) // 11 homozygote

{

mean[j] ;

cnt2[j] ;

mean2[j] ;

}

}

else

{

cnt2[j] ;

if ( ! *i2 ) // 00 homozygote

{

mean[j]--;

mean2[j]--;

}

}

}

else

{

if ( *i1 )

{

if ( *i2 ) // 11 homozygote

{

mean[j] ;

cnt1[j] ;

mean1[j] ;

}

}

else

{

cnt1[j] ;

if ( ! *i2 ) // 00 homozygote

{

mean[j]--;

mean1[j]--;

}

}

}

}

// Next individual

gperson ;

i1 ;

i2 ;

i ;

}

// Next SNP in set

ps ;

j ;

}

// Having iterated over all individuals, we can now calculate the mean

// values, perform mean-substitution of missing data, and calculate the

// second order terms

cout.precision(8);

for (int j=0; j<nss; j )

{

mean[j] /= (double)(cnt1[j] cnt2[j]);

mean1[j] /= (double)cnt1[j];

mean2[j] /= (double)cnt2[j];

}

// Element of pooled covariance matrix S is

// S[x][y] = ( X[i] - mean1[i] ) ( X[j] - mean1[j] )

pooled.resize(nss);

for (int j=0; j<nss; j )

pooled[j].resize(nss,0);

/////////////////////////////////////

// Iterate over pairs of SNPs in SET

// First SNP

vector<CSNP*>::iterator ps1 = pSNP.begin();

int j1=0;

while ( ps1 != pSNP.end() )

{

// Second SNP

vector<CSNP*>::iterator ps2 = ps1;

int j2=j1;

while ( ps2 != pSNP.end() )

{

///////////////////////////

// Iterate over individuals

vector<Individual*>::iterator gperson = sample.begin();

vector<bool>::iterator i1_1 = (*ps1)->one.begin();

vector<bool>::iterator i2_1 = (*ps1)->two.begin();

vector<bool>::iterator i1_2 = (*ps2)->one.begin();

vector<bool>::iterator i2_2 = (*ps2)->two.begin();

while ( gperson != sample.end() )

{

// Permuted self

Individual * pperson = (*gperson)->pperson;

// Set both values to sample mean

double v1 = mean[j1];

double v2 = mean[j2];

// First SNP

if ( *i1_1 )

{

if ( *i2_1 ) // 11 homozygote

{

v1 = 1;

}

}

else

{

if ( ! *i2_1 ) // 00 homozygote

{

v1 = -1;

}

else

v1 = 0; // 01 heterozygote

}

// Second SNP

if ( *i1_2 )

{

if ( *i2_2 ) // 11 homozygote

{

v2 = 1;

}

}

else

{

if ( ! *i2_2 ) // 00 homozygote

{

v2 = -1;

}

else

v2 = 0; // 01 heterozygote

}

// Contribution to covariance term

if (! pperson->missing)

{

if (pperson->aff) // affecteds

pooled[j1][j2] = ( v1 - mean2[j1] ) * ( v2 - mean2[j2] );

else // unaffecteds

pooled[j1][j2] = ( v1 - mean1[j1] ) * ( v2 - mean1[j2] );

}

// Next individual

gperson ;

i1_1 ;

i2_1 ;

i1_2 ;

i2_2 ;

}

// Next second SNP

ps2 ;

j2 ;

}

// Next first SNP

ps1 ;

j1 ;

}

// Make matrix symmetric

for (int i=0; i<nss; i )

for (int j=i; j<nss; j )

pooled[j][i] = pooled[i][j];

return;