Modified files for strain-specific residual variance option #269

jb621-star · 2023-03-20T14:12:36Z

These changes include the modifications to the math inside mvlmm, and attempts to add an option for the user to input their own strain-specific residual variance column within the files gemma_io.cpp and param.cpp. Thanks for your help!

… option similar to cvt

…eplaced in the math with the null Ve

Small formatting change to -residvar

Updated resid_var math to CalcQi

Updated most functions except MphEM and remaining P-val functions

Continued to make appropriate math changes to MphCalcSigma, UpdateV, MphCalcLogL, MphEM, MphCalcP, MphCalcBeta, and other functions.

Finished theoretical changes in other functions and replaced all "eps_eval" entries with "eval_vec, sigmasq".

Replaced any lingering mentions of "eps_eval" with "eval_vec, sigmasq".

Changed mentions of "eps_eval" with "sigmasq". Still need to figure out how to get eigenvectors from the same place "eval" is derived.

Changed mention of eps_eval as a vector to "gsl_matrix *sigmasq"

Changed the one mention of "eps_eval" to "eval_vec" and "sigmasq". May actually need to return to this file in order for the code to actually run.

eps_eval --> eval_vec, sigmasq

Changed eps_eval to sigmasq and corrected any identified regions where dimensions did not match.

Replaced input of CopyResid from eps_eval to sigmasq

aartaka

Here's what I've noticed. Most of the errors you're getting are due to mismatch between the arguments function requires and the arguments you provide, so ensure these are in sync. The rest I've tried to locate and comment on. Let's have a call (coordinate via Matrix) in case there's something I missed.

aartaka · 2024-10-20T04:25:48Z

src/param.h

  string file_anno; // Optional.
  string file_gxe;  // Optional.
  string file_cvt;  // Optional.
+  string file_residvar; // Optional


Notice you called it file_residvar here, but are using it as file_resid elsewhere. A typo?

src/gemma.cpp

aartaka · 2024-10-20T04:27:33Z

src/gemma.cpp

            // PLINK analysis
            if (cPar.file_gxe.empty()) {
-              cLmm.AnalyzePlink(U, eval, UtW, &UtY_col.vector, W,
+              cLmm.AnalyzePlink(U, eval, eval_vec, sigmasq, UtW, &UtY_col.vector, W,


You didn't define eval_vec anywhere earlier, that's why compiler is shouting. What is the intended use for eval_vec?

Also note that sigmasq is not defined either.

And then there's also required/provided argument mismatch. AnalyzePlink requires (U, eval, UtW, Uty, W, y, gwasnps), while you provide (U, eval, eval_vec, sigmasq, UtW, Uty, W, y, gwasnps). I guess AnalyzePlink needs to be modified to accept new arguments first.

aartaka · 2024-10-20T04:33:02Z

src/mvlmm.h


  string file_bfile;
  string file_geno;
+  //string file_residvar maybe?


Yes, if you need to use it (although you seems to only use residuals as sigmasq, so maybe no need for it.) And don't forget to add a line assigning it to CopyFromParam in case you add this field.

Changed "file_residvar" back to "file_resid"

Transferred all instances of "eval_vec" to "U".

Changed all instances of "_residvar" to "resid" to keep things consistent.

Changed all instances of "residvar" to "resid" to keep things consistent

Changed inputs to "CalcLmmVgVeBeta", "CalcLambda", and "AnalyzePlink" according to their use elsewhere (i.e. no more "eval_vec").

Changed residvar to resid for consistency.

Changed all mentions of residvar to resid for consistency

Changed all residvar to resid for consistency.

Changed "residvar" to "resid" and removed unnecessary mentions of "eval_vec" from other functions.

Redefined ReadFile_resid() to expect a gsl_matrix instead of <vector<vector<double>> for resid. This seems to be a consistent sticking point.

Changed ReadFile_resid so that gsl_matrix operations are being used for resid instead of vector<vector<double>> operations.

Threw out the following additional definition of ReadFile_resid: void ReadFile_resid(const string &file_resid, bool &error, gsl_matrix *sigmasq) { debug_msg("entered"); igzstream infile(file_resid.c_str(), igzstream::in); if (!infile) { cout << "error! fail to open the residual variance file: " << file_resid << endl; error = true; return; } size_t n_row = sigmasq->size1, n_col = sigmasq->size2, i_row = 0, i_col = 0; gsl_matrix_set_zero(sigmasq); //change this so that sigmasq = V_e somehow string line; char *ch_ptr; double d; while (getline(infile, line)) { if (i_row == n_row) { cout << "error! number of rows in the residual variance file is larger " << "than expected." << endl; error = true; } i_col = 0; ch_ptr = strtok((char *)line.c_str(), " ,\t"); while (ch_ptr != NULL) { if (i_col == n_col) { cout << "error! number of columns in the residual variance file " << "is larger than expected, for row = " << i_row << endl; error = true; } d = atof(ch_ptr); gsl_matrix_set(sigmasq, i_row, i_col, d); i_col++; ch_ptr = strtok(NULL, " ,\t"); } i_row++; } infile.close(); infile.clear(); return; }

Added additional ReadFile_resid call back, changing mentions of sigmasq back to resid.

Edited user documentation for "-resid" option.

Changed mention of "resid" to "sigmasq" within the function "CheckResid"

Created new version of "CheckResid" function that hopefully helps in debugging the empty "sigmasq" issue and commented out seemingly redundant "CheckResid" function initiation.

Restored ReadFile_resid to its original, because the other caused compile errors.

Copilot suggested a helper function for creating Sigma.

Suggested change to PARAM::CheckResid() according to Copilot.

Reinitialized Sigma in scope.

Modified some of the options here hoping that it would help me debug things a little better.

Updated PARAM::CheckResid to do a proper check in the case that n_resid is empty.

Added an additional debugging step for if the residual variance matrix is not initialized.

Corrected the CalculateSigma helper function for calculating the modified residual variance data structure.

Added various bug fixes with residual variance variables and checks to PARAM::CheckResid, PARAM::ReadFiles, and PARAM constructor.

Made sure to include Ve_null initiation.

Corrected subsequent function declarations of CalculateSigma.

Ve_null is never used, so removing from null initiation.

Removing Ve_null from initiation.

Setting resid to 1 if user doesn't supply residual variance file.

Attempted to address 'epsilon' warning by initializing it with all the other element-wise data structures.

Reordered PARAM::PARAM(void) parameters according to compiler warning messages.

Reordered PARAM terms according to compiler warnings.

Removed redeclaration of 'size_t PARAM::n_vc'

Adjustment to Matrix Copy Logic

Implemented fixes to CopyResid

Implemented warning to see if ni_test isn't set when CheckResid is run.

Called CheckResid AFTER ProcessCvtPhen() is called so that ni_test is already defined for default residal behavior.

jb621-star and others added 20 commits March 20, 2023 09:52

Modified files for strain-specific residual variance option

7c70ec5

Fixed 'cvt' variable and missing 0 < geno < 2 checks

ca93fb2

Fixed errors except for MphInitial errors

763a0b5

Fixed parameter names, attempted to integrate residual variance as an…

18e15e3

… option similar to cvt

More bug fixes

72f5a24

Added PARAM 'CheckResidvar' so that if no vector is supplied, it is r…

f76326b

…eplaced in the math with the null Ve

Update gemma.cpp

58e370f

Small formatting change to -residvar

Update mvlmm.cpp

924daa5

Updated resid_var math to CalcQi

Update mvlmm.cpp

db72327

Updated most functions except MphEM and remaining P-val functions

Update mvlmm.cpp

363ec00

Continued to make appropriate math changes to MphCalcSigma, UpdateV, MphCalcLogL, MphEM, MphCalcP, MphCalcBeta, and other functions.

Update mvlmm.cpp

c40699b

Finished theoretical changes in other functions and replaced all "eps_eval" entries with "eval_vec, sigmasq".

Update gemma.cpp

96d5f17

Replaced any lingering mentions of "eps_eval" with "eval_vec, sigmasq".

Update gemma.cpp

112d499

Update gemma_io.cpp

6f1d2cd

Changed mentions of "eps_eval" with "sigmasq". Still need to figure out how to get eigenvectors from the same place "eval" is derived.

Update gemma_io.h

057d27a

Changed mention of eps_eval as a vector to "gsl_matrix *sigmasq"

Update lmm.cpp

a8352e8

Changed the one mention of "eps_eval" to "eval_vec" and "sigmasq". May actually need to return to this file in order for the code to actually run.

Update mvlmm.h

f858b30

eps_eval --> eval_vec, sigmasq

Update param.cpp

587c181

Changed eps_eval to sigmasq and corrected any identified regions where dimensions did not match.

Update param.h

e85449d

Update param.h

ad54957

Replaced input of CopyResid from eps_eval to sigmasq

aartaka reviewed Oct 20, 2024

View reviewed changes

jb621-star added 9 commits October 22, 2024 12:13

Update param.h

b0da2e4

Changed "file_residvar" back to "file_resid"

Update gemma.cpp

9580e68

Transferred all instances of "eval_vec" to "U".

Update gemma_io.cpp

1b9c986

Changed all instances of "_residvar" to "resid" to keep things consistent.

Update gemma_io.h

e193600

Changed all instances of "residvar" to "resid" to keep things consistent

Update lmm.cpp

8409e49

Changed inputs to "CalcLmmVgVeBeta", "CalcLambda", and "AnalyzePlink" according to their use elsewhere (i.e. no more "eval_vec").

Update param.h

06b5f62

Changed residvar to resid for consistency.

Update param.cpp

be69a55

Changed all mentions of residvar to resid for consistency

Update gemma.cpp

2b25992

Changed all residvar to resid for consistency.

Update mvlmm.h

e1d92d0

Changed "residvar" to "resid" and removed unnecessary mentions of "eval_vec" from other functions.

jb621-star added 30 commits December 3, 2024 16:12

Update gemma_io.h

66a5a37

Redefined ReadFile_resid() to expect a gsl_matrix instead of <vector<vector<double>> for resid. This seems to be a consistent sticking point.

Update gemma_io.cpp

b6bfb79

Changed ReadFile_resid so that gsl_matrix operations are being used for resid instead of vector<vector<double>> operations.

Update gemma_io.cpp

ea678d7

Added additional ReadFile_resid call back, changing mentions of sigmasq back to resid.

Update gemma.cpp

40b31bc

Edited user documentation for "-resid" option.

Update param.cpp

bbb458b

Changed mention of "resid" to "sigmasq" within the function "CheckResid"

Update gemma_io.cpp

84cfda9

Created new version of "CheckResid" function that hopefully helps in debugging the empty "sigmasq" issue and commented out seemingly redundant "CheckResid" function initiation.

Update gemma_io.cpp

65dab2e

Restored ReadFile_resid to its original, because the other caused compile errors.

Update mvlmm.cpp

a96087b

Copilot suggested a helper function for creating Sigma.

Update param.cpp

d26bec6

Suggested change to PARAM::CheckResid() according to Copilot.

Update mvlmm.cpp

a8c75e0

Reinitialized Sigma in scope.

Update param.cpp

1aa998f

Modified some of the options here hoping that it would help me debug things a little better.

Update param.cpp

93641f1

Updated PARAM::CheckResid to do a proper check in the case that n_resid is empty.

Update gemma_io.cpp

1d49d03

Added an additional debugging step for if the residual variance matrix is not initialized.

Update mvlmm.cpp

528859d

Corrected the CalculateSigma helper function for calculating the modified residual variance data structure.

Update param.cpp

d8a1684

Added various bug fixes with residual variance variables and checks to PARAM::CheckResid, PARAM::ReadFiles, and PARAM constructor.

Update param.h

3c93381

Made sure to include Ve_null initiation.

Update mvlmm.cpp

79e04f1

Corrected subsequent function declarations of CalculateSigma.

Update param.cpp

a5e6474

Ve_null is never used, so removing from null initiation.

Update param.h

3840627

Removing Ve_null from initiation.

Update param.cpp

0839d6f

Setting resid to 1 if user doesn't supply residual variance file.

Update mvlmm.cpp

45ae373

Attempted to address 'epsilon' warning by initializing it with all the other element-wise data structures.

Update param.cpp

48ac836

Reordered PARAM::PARAM(void) parameters according to compiler warning messages.

Update param.h

d9daaba

Reordered PARAM terms according to compiler warnings.

Update param.h

9e5f10f

Removed redeclaration of 'size_t PARAM::n_vc'

Update param.cpp

1945cee

Update param.cpp

d08496d

Adjustment to Matrix Copy Logic

Update param.cpp

1592bce

Implemented fixes to CopyResid

Update param.cpp

a1e2445

Implemented warning to see if ni_test isn't set when CheckResid is run.

Update param.cpp

c19d26c

Called CheckResid AFTER ProcessCvtPhen() is called so that ni_test is already defined for default residal behavior.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Modified files for strain-specific residual variance option #269

Modified files for strain-specific residual variance option #269

Uh oh!

jb621-star commented Mar 20, 2023

Uh oh!

aartaka left a comment

Uh oh!

aartaka Oct 20, 2024

Uh oh!

Uh oh!

aartaka Oct 20, 2024

Uh oh!

aartaka Oct 20, 2024

Uh oh!

aartaka Oct 20, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Modified files for strain-specific residual variance option #269

Are you sure you want to change the base?

Modified files for strain-specific residual variance option #269

Uh oh!

Conversation

jb621-star commented Mar 20, 2023

Uh oh!

aartaka left a comment

Choose a reason for hiding this comment

Uh oh!

aartaka Oct 20, 2024

Choose a reason for hiding this comment

Uh oh!

Uh oh!

aartaka Oct 20, 2024

Choose a reason for hiding this comment

Uh oh!

aartaka Oct 20, 2024

Choose a reason for hiding this comment

Uh oh!

aartaka Oct 20, 2024

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants