[Genabel-commits] r1608 - pkg/VariABEL/src/VARlib
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 13 10:24:33 CET 2014
Author: maksim
Date: 2014-02-13 10:24:32 +0100 (Thu, 13 Feb 2014)
New Revision: 1608
Modified:
pkg/VariABEL/src/VARlib/supplementary_functions.h
pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp
Log:
Implementing 2df test. Almost done. Few bugs left.
Modified: pkg/VariABEL/src/VARlib/supplementary_functions.h
===================================================================
--- pkg/VariABEL/src/VARlib/supplementary_functions.h 2014-02-12 10:49:41 UTC (rev 1607)
+++ pkg/VariABEL/src/VARlib/supplementary_functions.h 2014-02-13 09:24:32 UTC (rev 1608)
@@ -74,8 +74,9 @@
{
public:
- my_small_vector(double * vector_, unsigned long number_)
+ my_small_vector(double * vector_, unsigned long number_, bool delocate_=true)
{
+ delocate = delocate_;
vector=vector_;
number=number_;
}
@@ -83,7 +84,7 @@
my_small_vector(const my_small_vector& p)
{
number = p.number;
-
+ delocate = true;
vector = new double[number];
for(int i=0 ; i<number ; i++)
{
@@ -94,13 +95,13 @@
~my_small_vector(void)
{
- if(vector != NULL) delete[] vector;
+ if(vector != NULL && delocate) delete[] vector;
}
double * vector;
long number; //amount of cells in vector
+ bool delocate;
-
};
//all info for a SNP is here
Modified: pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp
===================================================================
--- pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp 2014-02-12 10:49:41 UTC (rev 1607)
+++ pkg/VariABEL/src/VARlib/var_homogeneity_test_C.cpp 2014-02-13 09:24:32 UTC (rev 1608)
@@ -17,6 +17,7 @@
#include "var_homogeneity_tests.h"
+#include "supplementary_functions.h"
#include <iostream>
#include <iomanip>
#include <string>
@@ -110,11 +111,11 @@
}
// 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];
-// }
+ 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];
+ }
@@ -152,21 +153,15 @@
//___________________________________________
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];
-
-// for(int i=0; i<2*nids_nona ; i++)
-// {
-// }
-
-// for(int i=0; i<nids_nona ; i++)
-// {
-// }
-
-
-// for(int i=0 ; i<(*p_ * nids_nona) ; i++)
-// {
-// Rprintf("design_matrix[%d]=%f\n", i, design_matrix[i]);
-// }
+
+ if(*analys_type_ == 4)
+ {
+ for(long unsigned int i=0; i<nids_nona*step ; i++) design_matrix_copy[i+nids_nona] = design_matrix[i +(p-step)*nids_nona];
+ }
+ else
+ {
+ for(long unsigned int i=0; i<nids_nona ; i++) design_matrix_copy[i+nids_nona] = design_matrix[i +(p-1)*nids_nona];
+ }
linear_regression(
/*input variables:*/ trait, design_matrix, p_, &nids_nona,
@@ -194,35 +189,65 @@
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,
+ double *squared_residuals_ = new double[nids_nona];
+ for(int i=0; i<nids_nona ; i++) {squared_residuals_[i] = residuals[i]*residuals[i];}
+
+ my_small_vector squared_residuals(squared_residuals_, nids_nona);
+
+ linear_regression(
+ /*input variables:*/ squared_residuals_, design_matrix_copy, &p_disp, &nids_nona,
/*return variables:*/ betas_disp, se_disp, residuals,
/*auxiliary variables:*/ qty, jpvt, qraux, work, v, x_for_ch2inv);
//___________________________________________
+
+ my_small_vector squared_residuals_adusted_on_snp(residuals, nids_nona, false);
- 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];
- }
-
+ Rprintf("\n\n*analys_type_=%d\n", *analys_type_);
+ Rprintf("betas_disp[%d]=%f\n", 0, betas_disp[0]);
+ Rprintf("se_disp[%d]=%f\n", 0, se_disp[0]);
+
+
+
+ // Calculate significance of SNP through:
+ // SSRr <- sum(z^2) <- SSR for restricted model
+ // SSRu <- sum(z_residuals^2) <- SSR for unrestricted model
+ // q <- 1
+ // k <- 0
+ // Fstat <- ((SSRr-SSRu)/q)/(SSRu/(N-k-1))
+ double SSRr = my_var(&squared_residuals)*nids_nona;
+ double SSRu = my_var(&squared_residuals_adusted_on_snp)*nids_nona;
+ int q = step;
+
+ Rprintf("SSRr=%f, SSRu=%f", SSRr, SSRu);
+
+ static double Fstat;
+ Fstat = ((SSRr-SSRu)/q)/(SSRu/(nids_nona-1));
+
+// if(*analys_type_ == 4)
+// {
+ for(unsigned i=0 ; i<step ; i++)
+ {
+ Rprintf("betas_disp[%d]=%f\n", i+1, betas_disp[i+1]);
+ Rprintf("se_disp[%d]=%f\n", i+1, se_disp[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];
+// }
+
+ chi2[0] = Fstat;
+ df[0] = nids_nona;
+ Rprintf("END!!!!!!!!!!!");
}
return;
More information about the Genabel-commits
mailing list