diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/gemma_io.cpp | 4 | ||||
| -rw-r--r-- | src/mathfunc.cpp | 5 | ||||
| -rw-r--r-- | src/param.cpp | 1 |
3 files changed, 8 insertions, 2 deletions
diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp index 2380d45..cd2a5d1 100644 --- a/src/gemma_io.cpp +++ b/src/gemma_io.cpp @@ -2,7 +2,7 @@ Genome-wide Efficient Mixed Model Association (GEMMA) Copyright © 2011-2017, Xiang Zhou Copyright © 2017, Peter Carbonetto - Copyright © 2017, Pjotr Prins + Copyright © 2017-2020, Pjotr Prins This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -1494,7 +1494,7 @@ bool BimbamKin(const string file_geno, const set<string> ksnps, gsl_vector_add_constant(geno, -1.0 * geno_mean); - if (k_mode == 2 && geno_var != 0) { + if (k_mode == 2 && geno_var != 0) { // centering gsl_vector_scale(geno, 1.0 / sqrt(geno_var)); } // set the SNP column ns_test diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp index 614da14..9043e00 100644 --- a/src/mathfunc.cpp +++ b/src/mathfunc.cpp @@ -169,6 +169,7 @@ void CenterMatrix(gsl_matrix *G) { } // Center the matrix G. +// Only used in vc void CenterMatrix(gsl_matrix *G, const gsl_vector *w) { double d, wtw; gsl_vector *Gw = gsl_vector_safe_alloc(G->size1); @@ -192,6 +193,7 @@ void CenterMatrix(gsl_matrix *G, const gsl_vector *w) { } // Center the matrix G. +// Only used in vc void CenterMatrix(gsl_matrix *G, const gsl_matrix *W) { gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2); gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2); @@ -233,6 +235,7 @@ void CenterMatrix(gsl_matrix *G, const gsl_matrix *W) { } // "Standardize" the matrix G such that all diagonal elements = 1. +// (only used by vc) void StandardizeMatrix(gsl_matrix *G) { double d = 0.0; vector<double> vec_d; @@ -260,11 +263,13 @@ void StandardizeMatrix(gsl_matrix *G) { double ScaleMatrix(gsl_matrix *G) { double d = 0.0; + // Compute mean of diagonal for (size_t i = 0; i < G->size1; ++i) { d += gsl_matrix_get(G, i, i); } d /= (double)G->size1; + // Scale the matrix using the diagonal mean if (d != 0) { gsl_matrix_scale(G, 1.0 / d); } diff --git a/src/param.cpp b/src/param.cpp index 31b7382..cc4290c 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -1765,6 +1765,7 @@ void PARAM::CalcS(const map<string, double> &mapRS2wA, gsl_matrix_view Ksub = gsl_matrix_submatrix(K, 0, i * ni_test, ni_test, ni_test); CenterMatrix(&Ksub.matrix); + // Scale the matrix G such that the mean diagonal = 1. ScaleMatrix(&Ksub.matrix); gsl_matrix_view Asub = |