[Genabel-commits] r1607 - in pkg/VariABEL: R src/ITERlib src/VARlib
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 12 11:49:41 CET 2014
Author: maksim
Date: 2014-02-12 11:49:41 +0100 (Wed, 12 Feb 2014)
New Revision: 1607
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 2df test. In process... :-)
Modified: pkg/VariABEL/R/var_test_gwaa.R
===================================================================
--- pkg/VariABEL/R/var_test_gwaa.R 2014-02-11 08:29:43 UTC (rev 1606)
+++ pkg/VariABEL/R/var_test_gwaa.R 2014-02-12 10:49:41 UTC (rev 1607)
@@ -33,7 +33,7 @@
else if(testname == "likelihood") {testname <- 2}
else if(testname == "kolmogorov_smirnov") {testname <- 3}
else if(testname == "svlm") {testname <- 4}
-else {stop(paste(testname, "is unsupportive type of test. Only levene, and svlm are supported."))}
+else {stop(paste(testname, "is unsupportive type of test. Only levene and svlm are supported."))}
if(analysis_type == "AAvsABvsBB") {analysis_type <- 0}
@@ -41,13 +41,13 @@
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."))}
+else {stop(paste(analysis_type, "is unsupportive type of analysis. Only AAvsABvsBB, AAvsABandBB, ABvsAAandBB, BBvsAAandAB and 2df are supported."))}
MAR <- 2
OUT <- "R"
FUN <- "variance_homogeneity_test_C_wrapper"
STEP <- 1
-if(analysis_type == "2df")
+if(analysis_type == 4)
{
STEP <- 2
}
@@ -56,6 +56,7 @@
if(is(formula, "formula"))
{
+ print("Formula=formula")
design_matrix <- model.matrix(formula, data = phenodata)
# design_matrix <- design_matrix_
trat_name <- as.character(formula[[2]])
@@ -78,8 +79,7 @@
design_matrix_df_with_nas[na_bool,] <- design_matrix_df
design_matrix_df <- design_matrix_df_with_nas
}
-
- if(analysis_type == "2df")
+ if(analysis_type == 4)
{
design_matrix_df$snp_AB <- 0 #will be filled into iterator
design_matrix_df$snp_BB <- 0 #will be filled into iterator
@@ -96,13 +96,16 @@
# else if(is(formula, "numeric") || is(formula, "integer") || is(formula, "double"))
else if(is.character(formula))
{
+ print("Formula=char")
trait <- phenodata[,formula]
idnum <- length(trait)
- if(analysis_type == "2df")
+ print(c("analysis_type=", analysis_type))
+ if(analysis_type == 4)
{
design_matrix_df <- data.frame(X.Intercept=rep(1,idnum), snp_AB=rep(0,idnum), snp_BB=rep(0,idnum))
+ print("design_matrix_df 2df")
}
else
{
@@ -140,7 +143,14 @@
else if (class(genodata)=="databel")
{
gtNrow <- dim(genodata)[1]
- gtNcol <- dim(genodata)[2]
+ if(analysis_type == 4)
+ {
+ gtNcol <- dim(genodata)[2]/2
+ }
+ else
+ {
+ gtNcol <- dim(genodata)[2]
+ }
genodata <- genodata at data
if(!is.null(genodata_info))
@@ -224,6 +234,7 @@
#print(genodata_info_df)
#print(results_C)
+
results_df <- data.frame(genodata_info_df, results_C)
#print(results_df)
@@ -231,7 +242,7 @@
#print(cov_names)
-output_column_names <- c(colnames(genodata_info_df), "chisq", "df", "Intercept_effect", "Intercept_sd")
+output_column_names <- c(colnames(genodata_info_df), "chisq", "df", "Intercept_effect", "Intercept_se")
#print(output_column_names)
@@ -242,9 +253,23 @@
output_column_names <- c(output_column_names, paste(cov_names[i], "_se", sep=""))
}
-output_column_names <- c(output_column_names, "snp_eff_dispertion")
-output_column_names <- c(output_column_names, "snp_se_dispertion")
+
+if(analysis_type == 4)
+ {
+ output_column_names <- c(output_column_names, "snp_AB_eff_dispertion")
+ output_column_names <- c(output_column_names, "snp_AB_se_dispertion")
+
+ output_column_names <- c(output_column_names, "snp_BB_eff_dispertion")
+ output_column_names <- c(output_column_names, "snp_BB_se_dispertion")
+ }
+else
+ {
+ output_column_names <- c(output_column_names, "snp_eff_dispertion")
+ output_column_names <- c(output_column_names, "snp_se_dispertion")
+ }
+
+
colnames(results_df) <- output_column_names
if(testname != 4)
Modified: pkg/VariABEL/src/ITERlib/iterator.cpp
===================================================================
--- pkg/VariABEL/src/ITERlib/iterator.cpp 2014-02-11 08:29:43 UTC (rev 1606)
+++ pkg/VariABEL/src/ITERlib/iterator.cpp 2014-02-12 10:49:41 UTC (rev 1607)
@@ -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
@@ -337,6 +337,7 @@
**/
// Apply function of choosing
+ Rprintf("Run analysis\n");
pMethod(internal_data, nrow, step, out_data, ncol_multi, nrow_new, narg,
argList);
@@ -351,6 +352,7 @@
// = out_data[k];
//REAL(out)[(i_cs * ncol_multi + j) * nrow_new + k]
// = out_data[j*nrow_new+k];
+ Rprintf("Write ouput\n");
REAL(out)[j*ncol_cs+i_cs+k]
= out_data[j*nrow_new+k];
}
Modified: pkg/VariABEL/src/ITERlib/iterator_functions.cpp
===================================================================
--- pkg/VariABEL/src/ITERlib/iterator_functions.cpp 2014-02-11 08:29:43 UTC (rev 1606)
+++ pkg/VariABEL/src/ITERlib/iterator_functions.cpp 2014-02-12 10:49:41 UTC (rev 1607)
@@ -504,7 +504,7 @@
outdata[1] = double(*df);
int counter=2;
- for(int i=0 ; i<((*p)+1) ; i++)
+ for(int i=0 ; i<((*p)+step) ; i++)
{
outdata[counter] = betas[i];
counter++;
@@ -525,7 +525,7 @@
}
else
{
- outdataNcol = 4 + (*INTEGER(argList[3]))*2;
+ outdataNcol = 2 + step*2 + (*INTEGER(argList[3]))*2;
outdataNrow = 1;
}
Modified: pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp
===================================================================
--- pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp 2014-02-11 08:29:43 UTC (rev 1606)
+++ pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp 2014-02-12 10:49:41 UTC (rev 1607)
@@ -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, unsigned step)
+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 int step)
{
std::list<my_small_vector> trait_groups;
@@ -96,12 +96,26 @@
//fillfull the last covariates by the values from snp
- for(int i=nids*(p-1) ; i<nids*p ; i++)
- {
+
+ int start_snp;
+ int stop_snp = nids*p;
+
+ if(*analys_type_ == 4)
+ {
+ start_snp = nids*(p-2);
+ }
+ else
+ {
+ start_snp = nids*(p-1);
+ }
+
+// Rprintf("start_snp=%d, stop_snp=%d", start_snp, stop_snp);
+// for(int i=start_snp ; i<stop_snp; i++)
+// {
+// Rprintf("snp[%d]=%f, ", i-start_snp, snp[i-start_snp]);
+// design_matrix_copy[i]=snp[i-start_snp];
+// }
- design_matrix_copy[i]=snp[i-nids*(p-1)];
- }
-
@@ -149,12 +163,21 @@
// }
+// for(int i=0 ; i<(*p_ * nids_nona) ; i++)
+// {
+// Rprintf("design_matrix[%d]=%f\n", i, design_matrix[i]);
+// }
linear_regression(
/*input variables:*/ trait, design_matrix, p_, &nids_nona,
/*return variables:*/ betas, se, residuals,
/*auxiliary variables:*/ qty, jpvt, qraux, work, v, x_for_ch2inv);
+
+// for(int i=0 ; i<p ; i++)
+// {
+// Rprintf("betas[%d]=%f\n", i, betas[i]);
+// }
//___________________________________________
// double beta_snp = betas[1];
More information about the Genabel-commits
mailing list