[Genabel-commits] r1605 - in pkg/VariABEL: R src/ITERlib src/VARlib
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 11 08:55:37 CET 2014
Author: maksim
Date: 2014-02-11 08:55:37 +0100 (Tue, 11 Feb 2014)
New Revision: 1605
Modified:
pkg/VariABEL/R/var_test_gwaa.R
pkg/VariABEL/src/ITERlib/iterator.cpp
pkg/VariABEL/src/ITERlib/iterator_functions.cpp
pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp
Log:
Implementing of 2df test. Commit to save the current changes.
Modified: pkg/VariABEL/R/var_test_gwaa.R
===================================================================
--- pkg/VariABEL/R/var_test_gwaa.R 2014-02-10 22:29:48 UTC (rev 1604)
+++ pkg/VariABEL/R/var_test_gwaa.R 2014-02-11 07:55:37 UTC (rev 1605)
@@ -40,13 +40,18 @@
else if(analysis_type == "AAvsABandBB") {analysis_type <- 1}
else if(analysis_type == "ABvsAAandBB") {analysis_type <- 2}
else if(analysis_type == "BBvsAAandAB") {analysis_type <- 3}
+else if(analysis_type == "2df") {analysis_type <- 4}
else {stop(paste(analysis_type, "is unsupportive type of analysis. Only AAvsABvsBB, AAvsABandBB, ABvsAAandBB, and BBvsAAandAB are supported."))}
MAR <- 2
OUT <- "R"
FUN <- "variance_homogeneity_test_C_wrapper"
+STEP <- 1
+if(analysis_type == "2df")
+ {
+ STEP <- 2
+ }
-
idnum <- 0
if(is(formula, "formula"))
@@ -74,8 +79,17 @@
design_matrix_df <- design_matrix_df_with_nas
}
- design_matrix_df$snp <- 0 #will be filled into iterator
- idnum <- dim(design_matrix_df)[1]
+ if(analysis_type == "2df")
+ {
+ design_matrix_df$snp_AB <- 0 #will be filled into iterator
+ design_matrix_df$snp_BB <- 0 #will be filled into iterator
+ }
+ else
+ {
+ design_matrix_df$snp <- 0 #will be filled into iterator
+ }
+
+ idnum <- dim(design_matrix_df)[1]
}
@@ -84,7 +98,17 @@
{
trait <- phenodata[,formula]
idnum <- length(trait)
- design_matrix_df <- data.frame(X.Intercept=rep(1,idnum), snp=rep(0,idnum))
+
+
+ if(analysis_type == "2df")
+ {
+ design_matrix_df <- data.frame(X.Intercept=rep(1,idnum), snp_AB=rep(0,idnum), snp_BB=rep(0,idnum))
+ }
+ else
+ {
+ design_matrix_df <- data.frame(X.Intercept=rep(1,idnum), snp=rep(0,idnum))
+ }
+
}
#design_matrix_geno_means_df <- data.frame(X.Intercept=rep(1,idnum), snp=rep(0,idnum))
@@ -166,12 +190,14 @@
print("Start variance analysis...")
+
+
results_C <- .Call("iterator", genodata,
as.integer(gtNrow), as.integer(gtNcol),
as.character(FUN),
as.character(OUT),
- as.integer(MAR),
- as.integer(1),
+ as.integer(MAR), #margin
+ as.integer(STEP), #nrstep
as.integer(18), #iterator additional inputa parameters number
as.double(trait),
as.double(data.matrix(design_matrix_df)),
@@ -179,9 +205,9 @@
as.integer(p),
as.integer(analysis_type),
as.integer(testname),
- double(p+1),#betas
- double(p+1),#se
- double(1),#chi2
+ double(p+STEP),#betas
+ double(p+STEP),#se
+ double(STEP),#chi2
integer(1),#df
double(idnum),#residuals
double(idnum),#qty
Modified: pkg/VariABEL/src/ITERlib/iterator.cpp
===================================================================
--- pkg/VariABEL/src/ITERlib/iterator.cpp 2014-02-10 22:29:48 UTC (rev 1604)
+++ pkg/VariABEL/src/ITERlib/iterator.cpp 2014-02-11 07:55:37 UTC (rev 1605)
@@ -251,7 +251,7 @@
// Get the dimensions of the output the function of our choosing will be giving
pMethod(0, nrow, step, 0, ncol_multi, nrow_new, narg, argList);
- //Rprintf("ncol_multi,nrow_new=%d,%d\n",ncol_multi,nrow_new);
+// Rprintf("ncol_multi,nrow_new=%d,%d\n",ncol_multi,nrow_new);
// Allocate vector
// Start output SEXP for passing to R
// Even when the output is put into a filevector, we still return an (empty) SEXP
Modified: pkg/VariABEL/src/ITERlib/iterator_functions.cpp
===================================================================
--- pkg/VariABEL/src/ITERlib/iterator_functions.cpp 2014-02-10 22:29:48 UTC (rev 1604)
+++ pkg/VariABEL/src/ITERlib/iterator_functions.cpp 2014-02-11 07:55:37 UTC (rev 1605)
@@ -415,34 +415,13 @@
void variance_homogeneity_test_C_wrapper(double *indata, unsigned long int indataHeight,
- unsigned long int indataWidth,
+ unsigned long int step,
double *outdata, unsigned long int &outdataNcol,
unsigned long int &outdataNrow, unsigned int narg, SEXP *argList)
{
-// static double *design_matrix_new = new double[indataHeight*(*p)];
-// static double *design_matrix_copy_new = new double[indataHeight*(*p)];
-
-// SEXP design_matrix_new_SEXP;
-// SEXP design_matrix_copy_new_SEXP;
-// PROTECT(design_matrix_new_SEXP = allocVector(REALSXP, indataHeight*(*p)));
-// PROTECT(design_matrix_copy_new_SEXP = allocVector(REALSXP, indataHeight*(*p)));
-
-
-// double *design_matrix_new = REAL(design_matrix_new_SEXP);
-// double *design_matrix_copy_new = REAL(design_matrix_copy_new_SEXP);
-// UNPROTECT(2);
-
-
-// for(int i=0 ; i<indataHeight*(*p) ; i++)
-// {
-// design_matrix_new[i] = design_matrix[i];
-// design_matrix_copy_new[i] = design_matrix_copy[i];
-// }
-
-
if(indata)
{
double *trait_ = REAL(argList[0]);
@@ -517,7 +496,7 @@
variance_homogeneity_test_C(indata, trait, design_matrix_copy, design_matrix, p, &indataHeight, betas, se, chi2, df, residuals,
- analys_type, is_trait_na, testname, /*auxiliary variables:*/ qty, jpvt, qraux, work, v, x_for_ch2inv);
+ analys_type, is_trait_na, testname, /*auxiliary variables:*/ qty, jpvt, qraux, work, v, x_for_ch2inv, step);
// delete[] design_matrix_;
Modified: pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp
===================================================================
--- pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp 2014-02-10 22:29:48 UTC (rev 1604)
+++ pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp 2014-02-11 07:55:37 UTC (rev 1605)
@@ -66,7 +66,7 @@
//_________________________________________________________________________________________________
-void variance_homogeneity_test_C(double* snp, double *trait, double *design_matrix_copy, double *design_matrix, int *p_, long unsigned* nids_, double *betas, double *se, double * chi2, int * df, double * residuals, int *analys_type_, int * is_trait_na, int *testname, /*auxiliary variables:*/ double *qty, int *jpvt, double *qraux, double *work, double *v, double *x_for_ch2inv)
+void variance_homogeneity_test_C(double* snp, double *trait, double *design_matrix_copy, double *design_matrix, int *p_, long unsigned* nids_, double *betas, double *se, double * chi2, int * df, double * residuals, int *analys_type_, int * is_trait_na, int *testname, /*auxiliary variables:*/ double *qty, int *jpvt, double *qraux, double *work, double *v, double *x_for_ch2inv, unsigned step)
{
std::list<my_small_vector> trait_groups;
@@ -102,28 +102,9 @@
design_matrix_copy[i]=snp[i-nids*(p-1)];
}
-// for(int i=0 ; i<nids*p ; i++)
-// {
-// }
-
-// for(int i=0 ; i<nids ; i++)
-// {
-// }
-
-
-
- //drop out the NA values
- //___________________________________________
-// long unsigned int nids_nona=0;
-// for(int i=0 ; i<nids ; i++)
-// {
-// nids_nona += is_trait_na[i];
-// }
-// nids_nona = nids - nids_nona;
-
long unsigned int nids_nona=0;
for(int column=0 ; column<p ; column++)
{
@@ -156,7 +137,6 @@
//Get rid of major SNP effect
//___________________________________________
- int p_disp = 2;
for(long unsigned int i=0; i<nids_nona ; i++) design_matrix_copy[i] = design_matrix[i];
for(long unsigned int i=0; i<nids_nona ; i++) design_matrix_copy[i+nids_nona] = design_matrix[i +(p-1)*nids_nona];
@@ -187,8 +167,11 @@
//___________________________________________
- double betas_disp[2], se_disp[2];
+ //int p_disp = 2;
+ int p_disp = step + 1;
+ double betas_disp[step + 1], se_disp[step + 1];
+
for(int i=0; i<nids ; i++) {residuals[i] = residuals[i]*residuals[i];}
linear_regression(
/*input variables:*/ residuals, design_matrix_copy, &p_disp, &nids_nona,
@@ -198,14 +181,26 @@
//___________________________________________
-
-
- *chi2 = (betas_disp[1]/se_disp[1])*(betas_disp[1]/se_disp[1]);
- *df = 1;
- betas[p] = betas_disp[1];
- se[p] = se_disp[1];
+ if(*analys_type_ == 4)
+ {
+ for(unsigned i=0 ; i<2 ; i++)
+ {
+ chi2[i] = (betas_disp[i+1]/se_disp[i+1])*(betas_disp[i+1]/se_disp[i+1]);
+ df[i] = 1;
+ betas[p+i] = betas_disp[i+1];
+ se[p+i] = se_disp[i+1];
+ }
+ }
+ else
+ {
+ chi2[0] = (betas_disp[1]/se_disp[1])*(betas_disp[1]/se_disp[1]);
+ df[0] = 1;
+ betas[p] = betas_disp[1];
+ se[p] = se_disp[1];
+ }
- }
+
+ }
return;
} // end of function bartlett_test
@@ -213,170 +208,6 @@
-
-
-
-
} // end of extern "C"
-//extern "C" {
-//
-//
-//SEXP variance_homogeneity_test_C_old_data_type(SEXP set_, SEXP trait_, SEXP nids_, SEXP nsnps_, SEXP analys_type_, SEXP testname_)
-//{
-//char *set = (char*)(RAW(set_));
-//double *trait = REAL(trait_);
-//unsigned nids = INTEGER_VALUE(nids_);
-//int nsnps = INTEGER_VALUE(nsnps_);
-//const char * analys_type_char = CHARACTER_VALUE(analys_type_);
-//const char * testname_char = CHARACTER_VALUE(testname_);
-//
-//
-//
-//gtps_container Set(set, NULL, NULL, nids, nsnps); //creat object to facilitate working with set1
-//
-//int *snp = new int[nids];
-//double chi2;
-//double *chi2_vec;
-//int df;
-//int *df_vec;
-//
-//chi2_vec = new double[nsnps];
-//df_vec = new int[nsnps];
-//int* is_trait_na = new int[nids];
-//
-//
-////Fill is_trait_na. 1 means missing vallue, 0 - non missing
-////_________________________________________________________
-//for(unsigned id_counter=0 ; id_counter<nids ; id_counter++)
-// {
-// if(ISNAN(trait[id_counter]))
-// {
-// is_trait_na[id_counter]=1;
-// }
-// else
-// {
-// is_trait_na[id_counter]=0;
-// }
-// }
-////_________________________________________________________
-//
-//
-////_________________________________________________________
-//analysis_type_enum analys_type;
-//
-//if(strcmp(analys_type_char, "AAvsABvsBB") == 0)
-// {
-// analys_type = AAvsABvsBB;
-// }
-//else if(strcmp(analys_type_char, "AAvsABandBB") == 0)
-// {
-// analys_type = AAvsABandBB;
-// }
-//else if(strcmp(analys_type_char, "ABvsAAandBB") == 0)
-// {
-// analys_type = ABvsAAandBB;
-// }
-//else if(strcmp(analys_type_char, "BBvsAAandAB") == 0)
-// {
-// analys_type = BBvsAAandAB;
-// }
-//else
-// {
-// error("Unknown name of analysis type. Available are AAvsABvsBB, AAvsABandBB, ABvsAAandBB, BBvsAAandAB");
-// }
-////_________________________________________________________
-//
-//
-//
-////_________________________________________________________
-//testname_enum testname;
-//if(strcmp(testname_char, "bartlett") == 0)
-// {
-// testname = bartlett;
-// }
-//else if(strcmp(testname_char, "levene") == 0)
-// {
-// testname = levene;
-// }
-//else if(strcmp(testname_char, "likelihood") == 0)
-// {
-// testname = likelihood;
-// }
-//else if(strcmp(testname_char, "kolmogorov_smirnov") == 0)
-// {
-// testname = kolmogorov_smirnov;
-// }
-//else
-// {
-// error("Unknown name of test name. Available are bartlett, levene, likelihood, kolmogorov_smirnov");
-// }
-////_________________________________________________________
-//
-//unsigned step=10000;
-//
-//
-//for(int snp_position=1 ; snp_position<=nsnps ; snp_position++)
-// {
-// for(int id_counter=1 ; id_counter<=nids ; id_counter++) {snp[id_counter-1] = int(Set.get(id_counter, snp_position));}
-// variance_homogeneity_test_C(snp, trait, &nids, &chi2, &df, analys_type, is_trait_na, testname);
-// chi2_vec[snp_position-1] = chi2;
-// df_vec[snp_position-1] = df;
-//
-// if(snp_position % step == 0)
-// {
-// Rprintf("%d SNPs done\n", snp_position);
-// if(snp_position >= step*5) step *= 5;
-// }
-// }
-//
-////void variance_homogeneity_test_C(int* snp, double *trait, unsigned* nids, double * chi2, double * df, char* analys_type_, int * is_trait_na, char* testname)
-//
-//
-//
-//SEXP chi2_df_results_R_LIST;
-//PROTECT(chi2_df_results_R_LIST = NEW_LIST(2));
-//
-//SET_VECTOR_ELT(chi2_df_results_R_LIST, 0, NEW_NUMERIC(nsnps));
-//SET_VECTOR_ELT(chi2_df_results_R_LIST, 1, NEW_INTEGER(nsnps));
-//
-//for(unsigned i = 0; i < nsnps; i++)
-// {
-// NUMERIC_POINTER(VECTOR_ELT(chi2_df_results_R_LIST, 0))[i] = chi2_vec[i];
-// }
-//
-//for(unsigned i = 0; i < nsnps; i++)
-// {
-// INTEGER_POINTER(VECTOR_ELT(chi2_df_results_R_LIST, 1))[i] = df_vec[i];
-// }
-//
-//UNPROTECT(1);
-//
-////SEXP chi2_results_R;
-////PROTECT(chi2_results_R = allocVector(REALSXP, 2*nsnps));
-////double *chi2_results_R_double = REAL(chi2_results_R);
-////
-////for(unsigned i = 0; i < nsnps; i++)
-//// {
-//// chi2_results_R_double[i] = chi2_vec[i];
-//// }
-////
-////for(unsigned i = nsnps; i < 2*nsnps; i++)
-//// {
-//// chi2_results_R_double[i] = df_vec[i];
-//// }
-////UNPROTECT(1);
-//
-//
-//delete snp;
-//delete chi2_vec;
-//delete df_vec;
-//delete is_trait_na;
-//
-//return(chi2_df_results_R_LIST);
-//}
-//
-//
-//
-//} //end of extern "C"
More information about the Genabel-commits
mailing list