Browse code

Trying a third version of the E-GADGETS fitness score.

mnodzenski authored on 14/12/2022 15:08:05
Showing 2 changed files

... ...
@@ -2,7 +2,7 @@ Package: epistasisGA
2 2
 Type: Package
3 3
 Title: An R package to identify multi-snp effects in nuclear family studies 
4 4
     using the GADGETS method
5
-Version: 1.1.2
5
+Version: 1.1.3
6 6
 Authors@R: c(person("Michael", "Nodzenski", 
7 7
     email = "[email protected]",role = c("aut", "cre")), 
8 8
     person("Juno", "Krahn", role = "ctb"))
... ...
@@ -1173,7 +1173,7 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1173 1173
     if (!E_pd){
1174 1174
         
1175 1175
         arma::vec sum_dif_vecs(chrom_size, fill::ones);
1176
-        arma::vec beta_exposure_prob_disease(x_orig.n_cols, fill::zeros);
1176
+        arma::vec beta_prob_disease(x_orig.n_cols, fill::zeros);
1177 1177
         
1178 1178
         List res;
1179 1179
         res = List::create(Named("fitness_score") =  pow(10, -10),
... ...
@@ -1182,7 +1182,7 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1182 1182
                            Named("wald_stat") = pow(10, -10), 
1183 1183
                            Named("risk_set_alleles") = sum_diffs, 
1184 1184
                            Named("beta_exposure_prob_disease") = 
1185
-                               beta_exposure_prob_disease.t());
1185
+                               beta_prob_disease.t());
1186 1186
         return(res);
1187 1187
         
1188 1188
     } else {
... ...
@@ -1193,7 +1193,7 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1193 1193
         arma::vec prob_disease(mom_target.n_rows);
1194 1194
 
1195 1195
         //nominate risk alleles
1196
-        arma::mat pred_means = x_orig*beta_full;
1196
+        //arma::mat pred_means = x_orig*beta_full;
1197 1197
         for (unsigned int i = 0; i < mom_target.n_rows; i++){
1198 1198
             
1199 1199
             arma::vec risk_geno_probs(mom_target.n_cols);
... ...
@@ -1208,11 +1208,11 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1208 1208
                     
1209 1209
                     // zero out predicted mean difference vector for inconsistencies
1210 1210
                     // with nominated risk alleles
1211
-                    if (pred_means(i, j) <= 0){
1212
-                        
1213
-                        pred_means(i, j) = 0.0;
1214
-                        
1215
-                    }
1211
+                    // if (pred_means(i, j) <= 0){
1212
+                    //     
1213
+                    //     pred_means(i, j) = 0.0;
1214
+                    //     
1215
+                    // }
1216 1216
                     
1217 1217
                     if (mom_geno == 2.0 | dad_geno == 2.0){
1218 1218
                         
... ...
@@ -1236,11 +1236,11 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1236 1236
                     
1237 1237
                     // zero out predicted mean difference vector for inconsistencies
1238 1238
                     // with nominated risk alleles
1239
-                    if (pred_means(i, j) > 0){
1240
-                        
1241
-                        pred_means(i, j) = 0.0;
1242
-                        
1243
-                    }
1239
+                    // if (pred_means(i, j) > 0){
1240
+                    //     
1241
+                    //     pred_means(i, j) = 0.0;
1242
+                    //     
1243
+                    // }
1244 1244
                     
1245 1245
                     if (mom_geno == 0.0 | dad_geno == 0.0){
1246 1246
                         
... ...
@@ -1270,53 +1270,53 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1270 1270
         }
1271 1271
         
1272 1272
         // get the lengths of the predicted mean vecs
1273
-        arma::vec mean_vec_lengths = arma::sqrt(arma::sum(arma::square(pred_means), 1));
1273
+        // arma::vec mean_vec_lengths = arma::sqrt(arma::sum(arma::square(pred_means), 1));
1274 1274
     
1275 1275
         // make sure we have unique values (possible for all to be zero)    
1276
-        arma::vec unique_vec_lengths = arma::unique(mean_vec_lengths);
1277
-        
1278
-        // return small value if only one predicted mean 
1279
-        if (unique_vec_lengths.n_elem == 1){
1280
-            
1281
-            arma::vec sum_dif_vecs(chrom_size, fill::ones);
1282
-            arma::vec beta_exposure_prob_disease(x_orig.n_cols, fill::zeros);
1283
-            List res = List::create(Named("fitness_score") =  pow(10, -10),
1284
-                                    Named("sum_dif_vecs") = sum_dif_vecs.t(),
1285
-                                    Named("ht_trace") = pow(10, -10),
1286
-                                    Named("wald_stat") = pow(10, -10), 
1287
-                                    Named("risk_set_alleles") = sum_diffs, 
1288
-                                    Named("beta_exposure_prob_disease") = 
1289
-                                        beta_exposure_prob_disease.t());
1290
-            return(res);
1291
-            
1292
-            
1293
-        }
1276
+        // arma::vec unique_vec_lengths = arma::unique(mean_vec_lengths);
1277
+        // 
1278
+        // // return small value if only one predicted mean 
1279
+        // if (unique_vec_lengths.n_elem == 1){
1280
+        //     
1281
+        //     arma::vec sum_dif_vecs(chrom_size, fill::ones);
1282
+        //     arma::vec beta_exposure_prob_disease(x_orig.n_cols, fill::zeros);
1283
+        //     List res = List::create(Named("fitness_score") =  pow(10, -10),
1284
+        //                             Named("sum_dif_vecs") = sum_dif_vecs.t(),
1285
+        //                             Named("ht_trace") = pow(10, -10),
1286
+        //                             Named("wald_stat") = pow(10, -10), 
1287
+        //                             Named("risk_set_alleles") = sum_diffs, 
1288
+        //                             Named("beta_exposure_prob_disease") = 
1289
+        //                                 beta_exposure_prob_disease.t());
1290
+        //     return(res);
1291
+        //     
1292
+        //     
1293
+        // }
1294 1294
     
1295
-        arma::mat x_vec_lengths = join_rows(x0_orig, mean_vec_lengths);
1296
-        arma::vec beta_prob_disease = solve(x_vec_lengths, prob_disease, solve_opts::fast);
1297
-        arma::vec beta_exposure_prob_disease = solve(x_orig, prob_disease, solve_opts::fast);
1295
+        // arma::mat x_vec_lengths = join_rows(x0_orig, mean_vec_lengths);
1296
+        // arma::vec beta_prob_disease = solve(x_vec_lengths, prob_disease, solve_opts::fast);
1297
+        arma::vec beta_prob_disease = solve(x_orig, prob_disease, solve_opts::fast);
1298 1298
 
1299 1299
         // make sure association is positive 
1300
-        bool pos_assoc = beta_prob_disease(1) > 0;
1301
-        if (! pos_assoc){
1302
-            
1303
-            arma::vec sum_dif_vecs(chrom_size, fill::ones);
1304
-            List res = List::create(Named("fitness_score") =  pow(10, -10),
1305
-                                    Named("sum_dif_vecs") = sum_dif_vecs.t(),
1306
-                                    Named("ht_trace") = pow(10, -10),
1307
-                                    Named("wald_stat") = pow(10, -10), 
1308
-                                    Named("risk_set_alleles") = sum_diffs, 
1309
-                                    Named("beta_exposure_prob_disease") = 
1310
-                                        beta_exposure_prob_disease.t());
1311
-            
1312
-            return(res);
1313
-            
1314
-        }
1300
+        // bool pos_assoc = beta_prob_disease(1) > 0;
1301
+        // if (! pos_assoc){
1302
+        //     
1303
+        //     arma::vec sum_dif_vecs(chrom_size, fill::ones);
1304
+        //     List res = List::create(Named("fitness_score") =  pow(10, -10),
1305
+        //                             Named("sum_dif_vecs") = sum_dif_vecs.t(),
1306
+        //                             Named("ht_trace") = pow(10, -10),
1307
+        //                             Named("wald_stat") = pow(10, -10), 
1308
+        //                             Named("risk_set_alleles") = sum_diffs, 
1309
+        //                             Named("beta_exposure_prob_disease") = 
1310
+        //                                 beta_exposure_prob_disease.t());
1311
+        //     
1312
+        //     return(res);
1313
+        //     
1314
+        // }
1315 1315
         
1316
-        arma::colvec resid_prob_disease = prob_disease - x_vec_lengths*beta_prob_disease;
1316
+        arma::colvec resid_prob_disease = prob_disease - x_orig*beta_prob_disease;
1317 1317
         double sig2 = arma::as_scalar(arma::trans(resid_prob_disease)*resid_prob_disease/
1318
-                                      (x_vec_lengths.n_rows - x_vec_lengths.n_cols));
1319
-        arma::mat vcov_beta_prob_disease = sig2 * arma::pinv(arma::trans(x_vec_lengths)*x_vec_lengths);
1318
+                                      (x_orig.n_rows - x_orig.n_cols));
1319
+        arma::mat vcov_beta_prob_disease = sig2 * arma::pinv(arma::trans(x_orig)*x_orig);
1320 1320
         
1321 1321
         // make sure cov is positive definite and return small score if not
1322 1322
         bool vcov_beta_prob_disease_pd = vcov_beta_prob_disease.is_sympd();
... ...
@@ -1331,7 +1331,7 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1331 1331
                                     Named("wald_stat") = pow(10, -10), 
1332 1332
                                     Named("risk_set_alleles") = sum_diffs, 
1333 1333
                                     Named("beta_exposure_prob_disease") = 
1334
-                                        beta_exposure_prob_disease.t());
1334
+                                        beta_prob_disease.t());
1335 1335
             
1336 1336
             return(res);
1337 1337
             
... ...
@@ -1364,7 +1364,7 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1364 1364
                                Named("wald_stat") = wald_test, 
1365 1365
                                Named("risk_set_alleles") = sum_diffs, 
1366 1366
                                Named("beta_exposure_prob_disease") = 
1367
-                                   beta_exposure_prob_disease.t());
1367
+                                   beta_prob_disease.t());
1368 1368
             return(res);
1369 1369
             
1370 1370
         } else {
... ...
@@ -1391,7 +1391,7 @@ List GxE_fitness_score_mvlm(NumericMatrix case_genetic_data_,
1391 1391
                                Named("wald_stat") = wald_test, 
1392 1392
                                Named("risk_set_alleles") = sum_diffs, 
1393 1393
                                Named("beta_exposure_prob_disease") = 
1394
-                                   beta_exposure_prob_disease.t());
1394
+                                   beta_prob_disease.t());
1395 1395
             return(res);
1396 1396
             
1397 1397
         }