1 /*
2 * GCTA: a tool for Genome-wide Complex Trait Analysis
3 *
4 * Implementations of functions for bivariate REML analysis
5 *
6 * 2012 by Jian Yang <jian.yang@uq.edu.au>
7 *
8 * This file is distributed under the GNU General Public
9 * License, Version 2. Please see the file COPYING for more
10 * details
11 */
12
13 #include "gcta.h"
14
fit_bivar_reml(string grm_file,string phen_file,string qcovar_file,string covar_file,string keep_indi_file,string remove_indi_file,string sex_file,int mphen,int mphen2,double grm_cutoff,double adj_grm_fac,int dosage_compen,bool m_grm_flag,bool pred_rand_eff,bool est_fix_eff,int reml_mtd,int MaxIter,vector<double> reml_priors,vector<double> reml_priors_var,vector<int> drop,bool no_lrt,double prevalence,double prevalence2,bool no_constrain,bool ignore_Ce,vector<double> & fixed_rg_val,bool bivar_no_constrain)15 void gcta::fit_bivar_reml(string grm_file, string phen_file, string qcovar_file, string covar_file, string keep_indi_file, string remove_indi_file, string sex_file, int mphen, int mphen2, double grm_cutoff, double adj_grm_fac, int dosage_compen, bool m_grm_flag, bool pred_rand_eff, bool est_fix_eff, int reml_mtd, int MaxIter, vector<double> reml_priors, vector<double> reml_priors_var, vector<int> drop, bool no_lrt, double prevalence, double prevalence2, bool no_constrain, bool ignore_Ce, vector<double> &fixed_rg_val, bool bivar_no_constrain) {
16 _bivar_reml = true;
17 _bivar_no_constrain = bivar_no_constrain;
18 no_lrt = true;
19 _fixed_rg_val = fixed_rg_val;
20 _reml_mtd = reml_mtd;
21 _reml_max_iter = MaxIter;
22 int i = 0, j = 0, k = 0;
23 bool grm_flag = (!grm_file.empty());
24 bool qcovar_flag = (!qcovar_file.empty());
25 bool covar_flag = (!covar_file.empty());
26 if (m_grm_flag) grm_flag = false;
27
28 // Read data
29 stringstream errmsg;
30 int qcovar_num = 0, covar_num = 0;
31 vector<string> phen_ID, qcovar_ID, covar_ID, grm_id, grm_files;
32 vector< vector<string> > phen_buf, qcovar, covar; // save individuals by column
33
34 if (grm_flag) {
35 read_grm(grm_file, grm_id, true, false, !(adj_grm_fac > -1.0));
36 update_id_map_kp(grm_id, _id_map, _keep);
37 grm_files.push_back(grm_file);
38 }
39 else if (m_grm_flag) {
40 read_grm_filenames(grm_file, grm_files, false);
41 for (i = 0; i < grm_files.size(); i++) {
42 read_grm(grm_files[i], grm_id, false, true, !(adj_grm_fac > -1.0));
43 update_id_map_kp(grm_id, _id_map, _keep);
44 }
45 }
46 read_phen(phen_file, phen_ID, phen_buf, mphen, mphen2);
47 update_id_map_kp(phen_ID, _id_map, _keep);
48 if (qcovar_flag) {
49 qcovar_num = read_covar(qcovar_file, qcovar_ID, qcovar, true);
50 update_id_map_kp(qcovar_ID, _id_map, _keep);
51 }
52 if (covar_flag) {
53 covar_num = read_covar(covar_file, covar_ID, covar, false);
54 update_id_map_kp(covar_ID, _id_map, _keep);
55 }
56 if (!keep_indi_file.empty()) keep_indi(keep_indi_file);
57 if (!remove_indi_file.empty()) remove_indi(remove_indi_file);
58 if (grm_flag) {
59 if (grm_cutoff>-1.0) rm_cor_indi(grm_cutoff);
60 if (!sex_file.empty()) update_sex(sex_file);
61 if (adj_grm_fac>-1.0) adj_grm(adj_grm_fac);
62 if (dosage_compen>-1) dc(dosage_compen);
63 _grm_N.resize(0, 0);
64 }
65
66 vector<string> uni_id;
67 map<string, int> uni_id_map;
68 map<string, int>::iterator iter;
69 for (i = 0; i < _keep.size(); i++) {
70 uni_id.push_back(_fid[_keep[i]] + ":" + _pid[_keep[i]]);
71 uni_id_map.insert(pair<string, int>(_fid[_keep[i]] + ":" + _pid[_keep[i]], i));
72 }
73 _n = _keep.size();
74 if (_n < 1) throw ("Error: no individuals are in common in the input files.");
75 cout << _n << " individuals are in common in these files." << endl;
76
77 // construct model terms
78 int n1 = 0, n2 = 0;
79 vector<string> ystr1(_n), ystr2(_n);
80 mphen--;
81 mphen2--;
82 for (i = 0; i < phen_ID.size(); i++) {
83 iter = uni_id_map.find(phen_ID[i]);
84 if (iter == uni_id_map.end()) continue;
85 if (phen_buf[i][mphen] != "NA" && phen_buf[i][mphen] != "-9") n1++;
86 if (phen_buf[i][mphen2] != "NA" && phen_buf[i][mphen2] != "-9") n2++;
87 ystr1[iter->second] = phen_buf[i][mphen];
88 ystr2[iter->second] = phen_buf[i][mphen2];
89 }
90
91 _n = n1 + n2;
92 _y = eigenVector::Zero(_n);
93 vector<int> nms1, nms2;
94 int two_tr_comm = 0;
95 for (i = 0, j = 0, k = 0; i < _keep.size(); i++) {
96 bool tr1_miss = true;
97 if (ystr1[i] != "NA" && ystr1[i] != "-9") {
98 (_y.segment(0, n1))(j) = atof(ystr1[i].c_str());
99 nms1.push_back(i);
100 j++;
101 tr1_miss = false;
102 }
103 bool tr2_miss = true;
104 if (ystr2[i] != "NA" && ystr2[i] != "-9") {
105 (_y.segment(n1, n2))(k) = atof(ystr2[i].c_str());
106 nms2.push_back(i);
107 k++;
108 tr2_miss = false;
109 }
110 if (!tr1_miss && !tr2_miss) two_tr_comm++;
111 }
112 eigenVector y1_tmp = (_y.segment(0, n1)).array() - (_y.segment(0, n1)).mean();
113 _y_Ssq = y1_tmp.squaredNorm() / (n1 - 1.0);
114 if (!(fabs(_y_Ssq) < 1e30)) throw ("Error: the phenotypic variance for trait 1 is infinite. Please check the missing data in your phenotype file. Missing values should be represented by \"NA\" or \"-9\".");
115 eigenVector y2_tmp = (_y.segment(n1, n2)).array() - (_y.segment(n1, n2)).mean();
116 _y2_Ssq = y2_tmp.squaredNorm() / (n2 - 1.0);
117 if (!(fabs(_y2_Ssq) < 1e30)) throw ("Error: the phenotypic variance for trait 2 is infinite. Please check the missing data in your phenotype file. Missing values should be represented by \"NA\" or \"-9\".");
118 cout << nms1.size() << " non-missing phenotypes for trait #1 and " << nms2.size() << " for trait #2" << endl;
119 if (!ignore_Ce) {
120 if (two_tr_comm == 0) {
121 ignore_Ce = true;
122 cout << "Note: the residual covariance component is ignored because no individuals were measured for both traits." << endl;
123 } else if ((double) two_tr_comm / (double) _keep.size() < 0.1) {
124 ignore_Ce = true;
125 cout << "Note: the residual covariance component is ignored because < 10% of individuals were measured for both traits." << endl;
126 }
127 }
128
129 _ncase = 0.0;
130 _ncase2 = 0.0;
131 eigenVector y1 = _y.segment(0, n1), y2 = _y.segment(n1, n2);
132 _flag_CC = check_case_control(_ncase, y1);
133 if (_flag_CC) cout << "for trait #1" << endl;
134 else prevalence = -1.0;
135 _flag_CC2 = check_case_control(_ncase2, y2);
136 if (_flag_CC2) cout << "for trait #2" << endl;
137 else prevalence2 = -1.0;
138 //if(flag_CC2!=_flag_CC) throw("Error: for a bivariate analysis, the two traits should be both quantitative or both binary.");
139 if ((_flag_CC && prevalence<-1) || (_flag_CC2 && prevalence2<-1)) cout << "Note: we can specify the disease prevalence by the option --reml-bivar-prevalence so that GCTA can transform the variance explained to the underlying liability scale." << endl;
140
141 int pos = 0;
142 _r_indx.clear();
143 _bivar_pos.resize(3);
144 if (grm_flag) {
145 for (i = 0; i < 3 + 3 - ignore_Ce; i++) _r_indx.push_back(i);
146 _Asp.resize(_r_indx.size());
147 for (i = 0; i < _r_indx.size(); i++) (_Asp[i]).resize(_n, _n);
148 if (!no_lrt) drop_comp(drop);
149 _bivar_pos[0].push_back(pos);
150 for (j = 0; j < n1; j++) {
151 (_Asp[pos]).startVec(j);
152 for (i = 0; i < n1; i++) (_Asp[pos]).insertBack(i, j) = _grm(_keep[nms1[i]], _keep[nms1[j]]);
153 }
154 pos++;
155
156 _bivar_pos[1].push_back(pos);
157 for (j = 0; j < n2; j++) {
158 (_Asp[pos]).startVec(j + n1);
159 for (i = 0; i < n2; i++) (_Asp[pos]).insertBack(i + n1, j + n1) = _grm(_keep[nms2[i]], _keep[nms2[j]]);
160 }
161 pos++;
162
163 _bivar_pos[2].push_back(pos);
164 for (j = 0; j < n1; j++) {
165 (_Asp[pos]).startVec(j);
166 for (i = 0; i < n2; i++) (_Asp[pos]).insertBack(i + n1, j) = _grm(_keep[nms2[i]], _keep[nms1[j]]);
167 }
168 for (j = 0; j < n2; j++) {
169 (_Asp[pos]).startVec(j + n1);
170 for (i = 0; i < n1; i++) (_Asp[pos]).insertBack(i, j + n1) = _grm(_keep[nms1[i]], _keep[nms2[j]]);
171 }
172 pos++;
173
174 for (j = 0; j < pos; j++) (_Asp[j]).finalize();
175 _grm.resize(0, 0);
176 }
177 else if (m_grm_flag) {
178 if (!sex_file.empty()) update_sex(sex_file);
179 for (i = 0; i < 3 * grm_files.size() + 3 - ignore_Ce; i++) _r_indx.push_back(i);
180 _Asp.resize(_r_indx.size());
181 for (i = 0; i < _r_indx.size(); i++) (_Asp[i]).resize(_n, _n);
182 if (!no_lrt) drop_comp(drop);
183 string prev_file = grm_files[0];
184 vector<string> prev_grm_id(grm_id);
185 cout << "There are " << grm_files.size() << " GRM file names specified in the file [" + grm_file + "]." << endl;
186 vector<int> kp;
187 for (k = 0; k < grm_files.size(); k++) {
188 cout << "Reading the GRM from the " << k + 1 << "th file ..." << endl;
189 read_grm(grm_files[k], grm_id, true, false, !(adj_grm_fac > -1.0));
190 if (adj_grm_fac>-1.0) adj_grm(adj_grm_fac);
191 if (dosage_compen>-1) dc(dosage_compen);
192 StrFunc::match(uni_id, grm_id, kp);
193 int pos0 = pos;
194
195 _bivar_pos[0].push_back(pos);
196 for (j = 0; j < n1; j++) {
197 (_Asp[pos]).startVec(j);
198 for (i = 0; i < n1; i++) (_Asp[pos]).insertBack(i, j) = _grm(kp[nms1[i]], _keep[nms1[j]]);
199 }
200 pos++;
201
202 _bivar_pos[1].push_back(pos);
203 for (j = 0; j < n2; j++) {
204 (_Asp[pos]).startVec(j + n1);
205 for (i = 0; i < n2; i++) (_Asp[pos]).insertBack(i + n1, j + n1) = _grm(kp[nms2[i]], _keep[nms2[j]]);
206 }
207 pos++;
208
209 _bivar_pos[2].push_back(pos);
210 for (j = 0; j < n1; j++) {
211 (_Asp[pos]).startVec(j);
212 for (i = 0; i < n2; i++) (_Asp[pos]).insertBack(i + n1, j) = _grm(kp[nms2[i]], _keep[nms1[j]]);
213 }
214 for (j = 0; j < n2; j++) {
215 (_Asp[pos]).startVec(j + n1);
216 for (i = 0; i < n1; i++) (_Asp[pos]).insertBack(i, j + n1) = _grm(kp[nms1[i]], _keep[nms2[j]]);
217 }
218 pos++;
219
220 for (j = pos0; j < pos; j++) (_Asp[j]).finalize();
221 prev_file = grm_files[k];
222 prev_grm_id = grm_id;
223 }
224 _grm_N.resize(0, 0);
225 _grm.resize(0, 0);
226 }
227
228 _bivar_pos[0].push_back(pos);
229 for (i = 0; i < n1; i++) {
230 (_Asp[pos]).startVec(i);
231 (_Asp[pos]).insertBack(i, i) = 1.0;
232 }
233 (_Asp[pos]).finalize();
234 pos++;
235
236 _bivar_pos[1].push_back(pos);
237 for (i = 0; i < n2; i++) {
238 (_Asp[pos]).startVec(i + n1);
239 (_Asp[pos]).insertBack(i + n1, i + n1) = 1.0;
240 }
241 (_Asp[pos]).finalize();
242 pos++;
243
244 if (!ignore_Ce) {
245 _bivar_pos[2].push_back(pos);
246 for (j = 0; j < n1; j++) {
247 (_Asp[pos]).startVec(j);
248 for (i = 0; i < n2; i++) {
249 if (nms2[i] == nms1[j]) (_Asp[pos]).insertBack(i + n1, j) = 1;
250 }
251 }
252 for (j = 0; j < n2; j++) {
253 (_Asp[pos]).startVec(j + n1);
254 for (i = 0; i < n1; i++) {
255 if (nms1[i] == nms2[j]) (_Asp[pos]).insertBack(i, j + n1) = 1;
256 }
257 }
258 pos++;
259 }
260
261 // construct X matrix
262 vector<eigenMatrix> E_float;
263 eigenMatrix qE_float;
264 construct_X(_keep.size(), uni_id_map, qcovar_flag, qcovar_num, qcovar_ID, qcovar, covar_flag, covar_num, covar_ID, covar, E_float, qE_float);
265 eigenMatrix X(_X);
266 _X = eigenMatrix::Zero(_n, _X_c * 2);
267 for (i = 0; i < n1; i++) (_X.block(0, 0, n1, _X_c)).row(i) = X.row(nms1[i]);
268 for (i = 0; i < n2; i++) (_X.block(n1, _X_c, n2, _X_c)).row(i) = X.row(nms2[i]);
269 _X_c *= 2;
270
271 // names of variance component
272 for (i = 0; i < grm_files.size(); i++) {
273 stringstream strstrm;
274 if (grm_files.size() == 1) strstrm << "";
275 else strstrm << i + 1;
276 _var_name.push_back("V(G" + strstrm.str() + ")_tr1");
277 _var_name.push_back("V(G" + strstrm.str() + ")_tr2");
278 _var_name.push_back("C(G" + strstrm.str() + ")_tr12");
279 _hsq_name.push_back("V(G" + strstrm.str() + ")/Vp_tr1");
280 _hsq_name.push_back("V(G" + strstrm.str() + ")/Vp_tr2");
281 }
282 _var_name.push_back("V(e)_tr1");
283 _var_name.push_back("V(e)_tr2");
284 if (!ignore_Ce) _var_name.push_back("C(e)_tr12");
285 _ignore_Ce = ignore_Ce;
286
287 // run REML algorithm
288 reml(pred_rand_eff, est_fix_eff, reml_priors, reml_priors_var, prevalence, prevalence2, no_constrain, no_lrt, false);
289 }
290
calcu_Vi_bivar(eigenMatrix & Vi,eigenVector & prev_varcmp,double & logdet,int & iter)291 bool gcta::calcu_Vi_bivar(eigenMatrix &Vi, eigenVector &prev_varcmp, double &logdet, int &iter) {
292 int i = 0, n = Vi.cols();
293 double d_buf = 0.0;
294 logdet = 0.0;
295 string errmsg = "\nError: the V (variance-covariance) matrix is not invertible.";
296
297 Vi = eigenMatrix::Zero(_n, _n);
298 for (i = 0; i < _r_indx.size(); i++) Vi += (_Asp[_r_indx[i]]) * prev_varcmp[i];
299
300 if (_V_inv_mtd == 0) {
301 if (!comput_inverse_logdet_LDLT_mkl(Vi, logdet)) {
302 //cout<<"Note: the variance-covaraince matrix V is non-positive definite. Switching to Cholesky to LU decomposition approach."<<endl;
303 _V_inv_mtd = 1;
304 }
305 }
306 if (_V_inv_mtd == 1) {
307 if (!comput_inverse_logdet_LU_mkl(Vi, logdet)) throw ("Error: the variance-covaraince matrix V is not invertible.");
308 }
309
310 /* if(!comput_inverse_logdet_LDLT_mkl(Vi, logdet)){
311 if(_reml_have_bend_A) throw(errmsg);
312 cout<<"Warning: the variance-covaraince matrix V is negative-definite."<<endl;
313 bend_A();
314 _reml_have_bend_A=true;
315 iter=-1;
316 cout<<"Restarting iterations ..."<<endl;
317 return false;
318 }
319 */
320
321 }
322
constrain_rg(eigenVector & varcmp)323 void gcta::constrain_rg(eigenVector &varcmp) {
324 static int count = 0;
325 int v_pos = 0, c_pos = 0, d = _bivar_pos[0].size() - _bivar_pos[2].size();
326 if (_ignore_Ce) d--;
327 eigenMatrix G(2, 2);
328
329 for (c_pos = _bivar_pos[2].size() - 1; c_pos >= 0; c_pos--) {
330 v_pos = c_pos + d;
331 G(0, 0) = varcmp[_bivar_pos[0][v_pos]];
332 G(1, 1) = varcmp[_bivar_pos[1][v_pos]];
333 G(0, 1) = G(1, 0) = varcmp[_bivar_pos[2][c_pos]];
334
335 SelfAdjointEigenSolver<eigenMatrix> eigensolver(G);
336 eigenVector eval = eigensolver.eigenvalues();
337 if (eval.minCoeff() <= 0.0) {
338 if (count == 0) {
339 cout << "Note: to constrain the correlation being from -1 to 1, a genetic (or residual) variance-covariance matrix is bended to be positive definite." << endl;
340 count++;
341 }
342 bending_eigenval(eval);
343 G = eigensolver.eigenvectors() * eigenDiagMat(eval) * eigensolver.eigenvectors().transpose();
344 varcmp[_bivar_pos[0][v_pos]] = G(0, 0);
345 varcmp[_bivar_pos[1][v_pos]] = G(1, 1);
346 varcmp[_bivar_pos[2][c_pos]] = G(0, 1);
347 }
348 }
349 }
350
calcu_rg(eigenVector & varcmp,eigenMatrix & Hi,eigenVector & rg,eigenVector & rg_var,vector<string> & rg_name)351 void gcta::calcu_rg(eigenVector &varcmp, eigenMatrix &Hi, eigenVector &rg, eigenVector &rg_var, vector<string> &rg_name) {
352 int i = 0, j = 0;
353 double V1 = 0, V2 = 0, C = 0, VarV1 = 0, VarV2 = 0, VarC = 0.0, CovV1V2 = 0.0, CovV1C = 0.0, CovV2C = 0.0;
354
355 rg = eigenVector::Zero(_bivar_pos[0].size() - 1);
356 rg_var = rg;
357 for (i = 0; i < _bivar_pos[0].size() - 1; i++) {
358 V1 = varcmp[_bivar_pos[0][i]];
359 V2 = varcmp[_bivar_pos[1][i]];
360 C = varcmp[_bivar_pos[2][i]];
361 VarV1 = Hi(_bivar_pos[0][i], _bivar_pos[0][i]);
362 VarV2 = Hi(_bivar_pos[1][i], _bivar_pos[1][i]);
363 VarC = Hi(_bivar_pos[2][i], _bivar_pos[2][i]);
364 CovV1V2 = Hi(_bivar_pos[0][i], _bivar_pos[1][i]);
365 CovV1C = Hi(_bivar_pos[0][i], _bivar_pos[2][i]);
366 CovV2C = Hi(_bivar_pos[1][i], _bivar_pos[2][i]);
367 if (V1 * V2 > 0) {
368 rg[i] = sqrt(V1 * V2);
369 if (rg[i] > 0) rg[i] = C / rg[i];
370 rg_var[i] = rg[i] * rg[i]*(VarV1 / (4 * V1 * V1) + VarV2 / (4 * V2 * V2) + VarC / (C * C) + CovV1V2 / (2 * V1 * V2) - CovV1C / (V1 * C) - CovV2C / (V2 * C));
371 }
372
373 if (_bivar_pos[0].size() == 2) rg_name.push_back("rG");
374 else {
375 stringstream strstrm;
376 strstrm << "rG" << i + 1;
377 rg_name.push_back(strstrm.str());
378 }
379 }
380 }
381
lgL_fix_rg(eigenVector & prev_varcmp,bool no_constrain)382 double gcta::lgL_fix_rg(eigenVector &prev_varcmp, bool no_constrain) {
383 int i = 0, j = 0;
384 if (_fixed_rg_val.size() > _bivar_pos[0].size() - 1) {
385 vector<double> rg_val_buf(_fixed_rg_val);
386 _fixed_rg_val.clear();
387 for (i = 0; i < _bivar_pos[0].size() - 1; i++) _fixed_rg_val.push_back(rg_val_buf[i]);
388 }
389
390 cout << "\nCalculating the logLikelihood for the model with the genetic correlation" << (_fixed_rg_val.size() > 1 ? "s" : "") << " being fixed at ";
391 for (int i = 0; i < _fixed_rg_val.size() - 1; i++) cout << _fixed_rg_val[i] << "\t";
392 cout << _fixed_rg_val[_fixed_rg_val.size() - 1] << endl;
393
394 vector<int> r_indx_buf(_r_indx);
395 _bivar_pos_prev = _bivar_pos;
396 for (i = _fixed_rg_val.size() - 1; i >= 0; i--) {
397 for (j = i + 1; j < _bivar_pos_prev[0].size(); j++) _bivar_pos[0][j]--;
398 for (j = i + 1; j < _bivar_pos_prev[1].size(); j++) _bivar_pos[1][j]--;
399 for (j = i + 1; j < _bivar_pos_prev[2].size(); j++) _bivar_pos[2][j]--;
400 _r_indx.erase(_r_indx.begin() + _bivar_pos[2][i]);
401 _bivar_pos[2].erase(_bivar_pos[2].begin() + i);
402 }
403
404 _Asp_prev = _Asp;
405 eigenMatrix Vi_X(_n, _X_c), Xt_Vi_X_i(_X_c, _X_c), Hi(_r_indx.size(), _r_indx.size());
406 eigenVector Py(_n);
407 eigenVector varcmp(_r_indx.size());
408 for (i = 0; i < _r_indx.size(); i++) {
409 varcmp[i] = fabs(prev_varcmp[_r_indx[i]]);
410 if (varcmp[i] < 1.0e-30) varcmp[i] = 0.1;
411 }
412 double lgL = reml_iteration(Vi_X, Xt_Vi_X_i, Hi, Py, varcmp, false, no_constrain, true);
413 _r_indx = r_indx_buf;
414 _bivar_pos = _bivar_pos_prev;
415
416 return lgL;
417 }
418
update_A(eigenVector & prev_varcmp)419 void gcta::update_A(eigenVector &prev_varcmp) {
420 int i = 0;
421 double g1 = 0.0, g2 = 0.0;
422 for (i = 0; i < _fixed_rg_val.size(); i++) {
423 g1 = prev_varcmp[_bivar_pos[0][i]];
424 g2 = prev_varcmp[_bivar_pos[1][i]];
425 _Asp[_bivar_pos_prev[0][i]] = _Asp_prev[_bivar_pos_prev[0][i]] + _Asp_prev[_bivar_pos_prev[2][i]]*(0.5 * _fixed_rg_val[i] * sqrt(g2 / g1));
426 _Asp[_bivar_pos_prev[1][i]] = _Asp_prev[_bivar_pos_prev[1][i]] + _Asp_prev[_bivar_pos_prev[2][i]]*(0.5 * _fixed_rg_val[i] * sqrt(g1 / g2));
427 }
428
429 }
430
431 /*
432 void gcta::init_rg(eigenVector &varcmp)
433 {
434 int i=0;
435 double d_buf=0.0, ratio=0.0, delta_0=0.0, delta_1=0.0;
436 eigenVector varcmp_buf=varcmp;
437
438 for(i=0; i<_fixed_rg_val.size(); i++){
439 d_buf=varcmp[_bivar_pos[0][i]]*varcmp[_bivar_pos[1][i]];
440 if(d_buf>0.0){
441 ratio=sqrt(abs(_fixed_rg_val[i]/((varcmp[_bivar_pos[2][i]]/sqrt(d_buf)))));
442 if(d_buf*_fixed_rg_val[i]>0.0) varcmp[_bivar_pos[2][i]]*=ratio;
443 else varcmp[_bivar_pos[2][i]]*=(-1.0*ratio);
444 varcmp[_bivar_pos[0][i]]/=ratio;
445 varcmp[_bivar_pos[1][i]]/=ratio;
446 delta_0=varcmp_buf[_bivar_pos[0][i]]-varcmp[_bivar_pos[0][i]];
447 }
448 else throw("Error: unable to calcuate the genetic correlation because one of the genetic variance components is negative.");
449 }
450 }
451
452 void gcta::fix_rg(eigenVector &varcmp, int pos)
453 {
454 int i=0;
455 double d_buf=0.0;
456 string errmsg="Error: unable to calcuate the genetic correlation because one of the genetic variance components is zero.";
457
458 for(i=0; i<_fixed_rg_val.size(); i++){
459 if(pos==2){
460 d_buf=varcmp[_bivar_pos[0][i]]*varcmp[_bivar_pos[1][i]];
461 if(CommFunc::FloatNotEqual(d_buf, 0.0)) varcmp[_bivar_pos[2][i]]=_fixed_rg_val[i]*sqrt(d_buf);
462 else throw(errmsg);
463 }
464 else if(pos==1){
465 if(CommFunc::FloatNotEqual(varcmp[_bivar_pos[1][i]], 0.0)) varcmp[_bivar_pos[0][i]]=varcmp[_bivar_pos[2][i]]*varcmp[_bivar_pos[2][i]]/(varcmp[_bivar_pos[1][i]]*_fixed_rg_val[i]*_fixed_rg_val[i]);
466 else throw(errmsg);
467 }
468 else if(pos==0){
469 if(CommFunc::FloatNotEqual(varcmp[_bivar_pos[0][i]], 0.0)) varcmp[_bivar_pos[1][i]]=varcmp[_bivar_pos[2][i]]*varcmp[_bivar_pos[2][i]]/(varcmp[_bivar_pos[0][i]]*_fixed_rg_val[i]*_fixed_rg_val[i]);
470 else throw(errmsg);
471 }
472 }
473 }
474 */