[Genabel-commits] r910 - pkg/ProbABEL/src

[noreply at r-forge.r-project.org]

Author: lckarssen
Date: 2012-05-21 22:06:03 +0200 (Mon, 21 May 2012)
New Revision: 910

Modified:
   pkg/ProbABEL/src/data.h
   pkg/ProbABEL/src/reg1.h
Log:
ProbABEL:
"Forward port" of the fix from rev. 902 in the ProbABEL-pacoxph branch where the behaviour for pacoxph was fixed for .prob files (ngpreds=2). 

WARNING: NOT WELL TESTED YET BECAUSE pacoxph STILL RUNS VERY SLOW SINCE THE ADDITION OF FILEVECTOR STUFF.



Modified: pkg/ProbABEL/src/data.h
===================================================================
--- pkg/ProbABEL/src/data.h	2012-05-21 19:29:21 UTC (rev 909)
+++ pkg/ProbABEL/src/data.h	2012-05-21 20:06:03 UTC (rev 910)
@@ -72,8 +72,8 @@
     int noutcomes;
     int ncov;
     unsigned short int * allmeasured;
-    mematrix<double> X;
-    mematrix<double> Y;
+    mematrix<double> X;		/* Will contain the values of the covariate(s) */
+    mematrix<double> Y;		/* Will contain the values of the outcome(s) */
     std::string * idnames;
     std::string model;
     std::string * model_terms;
@@ -98,7 +98,7 @@
 	    {
 
 		nphenocols++;
-		//				std::cout << tmp << " " << nphenocols << " ";
+		//   std::cout << tmp << " " << nphenocols << " ";
 	    }
 	    while (myfile.getline(line,BFS)) {
 		int tmplins = 0;
@@ -531,11 +531,14 @@
 	for (int j=0;j<ngpreds;j++)
 	{
 	    float snpdata[nids];
-	    for (int i=0;i<nids;i++) masked_data[i]=0;
-	    gend.get_var(snpnum*ngpreds+j,snpdata);
-	    for (int i=0;i<nids;i++) {
-		X.put(snpdata[i],i,(ncov+1-j-1));
-		if (isnan(snpdata[i])) masked_data[i]=1;
+	    for (int i=0; i<nids; i++) masked_data[i]=0;
+
+	    gend.get_var(snpnum*ngpreds+j, snpdata);
+
+	    for (int i=0; i<nids; i++)
+	    {
+		X.put(snpdata[i], i, (ncov-j));
+		if (isnan(snpdata[i])) masked_data[i] = 1;
 	    }
 	}
     }
@@ -672,7 +675,7 @@
 	    sstat[i] = int((phed.Y).get(i,1));
 	    if (sstat[i] != 1 && sstat[i]!=0)
 	    {
-		std::cerr << "coxph_data: status not 0/1 (right order: id, fuptime, status ...)"
+		std::cerr << "coxph_data: status not 0/1 (correct order: id, fuptime, status ...)"
 			  << endl;
 		exit(1);
 	    }
@@ -745,7 +748,12 @@
 
     void update_snp(gendata &gend, int snpnum)
     {
-	// note this sorts by "order"!!!
+	/** note this sorts by "order"!!!
+	 * Here we deal with transposed X, hence last two arguments are swapped
+	 * compared to the other 'update_snp'
+	 * Also, the starting column-1 is not necessary for cox X therefore
+	 * 'ncov-j' changes to 'ncov-j-1'
+	 **/
 
 	for (int j=0; j<ngpreds; j++)
 	{
@@ -756,14 +764,11 @@
 	    gend.get_var(snpnum*ngpreds+j, snpdata);
 
 	    for (int i=0; i<nids; i++) {
-		X.put(snpdata[i], (ncov-ngpreds+j), order[i]);
+		X.put(snpdata[i], (ncov-j-1), order[i]);
 		if ( isnan(snpdata[i]) )
 		    masked_data[order[i]] = 1;
 	    }
 	}
-	//		for (int i=0;i<nids;i++)
-	//			for (int j=0;j<ngpreds;j++)
-	//				X.put((gend.G).get(i,(snpnum*ngpreds+j)),(ncov-ngpreds+j),order[i]);
     }
 
     ~coxph_data()

Modified: pkg/ProbABEL/src/reg1.h
===================================================================
--- pkg/ProbABEL/src/reg1.h	2012-05-21 19:29:21 UTC (rev 909)
+++ pkg/ProbABEL/src/reg1.h	2012-05-21 20:06:03 UTC (rev 910)
@@ -173,21 +173,21 @@
     {
 	nX.reinit(X.nrow,(X.ncol-1));
     }
-    int c1 = X.ncol-2;
-    int c2 = X.ncol-1;
+    int c1 = X.ncol-2; 		/* column with Prob(A1A2) */
+    int c2 = X.ncol-1;		/* column with Prob(A1A1). Note the order is swapped cf the file! */
     for (int i=0; i<X.nrow; i++)
 	for (int j=0; j<(X.ncol-2); j++)
 	    nX[i * nX.ncol + j] = X[i * X.ncol + j];
 
     for (int i=0; i<nX.nrow; i++)
     {
-	if (model==1)
+	if (model==1)		/* Additive */
 	    nX[i*nX.ncol+c1] = X[i*X.ncol+c1] + 2.*X[i*X.ncol+c2];
-	else if (model==2)
+	else if (model==2)	/* Dominant */
 	    nX[i*nX.ncol+c1] = X[i*X.ncol+c1] + X[i*X.ncol+c2];
-	else if (model==3)
+	else if (model==3)	/* Recessive */
 	    nX[i*nX.ncol+c1] = X[i*X.ncol+c2];
-	else if (model==4)
+	else if (model==4)	/* over-dominant */
 	    nX[i*nX.ncol+c1] = X[i*X.ncol+c1];
 	if(interaction != 0)
 	    nX[i*nX.ncol+c2] = X[i*nX.ncol+interaction] * nX[i*nX.ncol+c1]; //Maksim: interaction with SNP
@@ -920,10 +920,9 @@
 		  double tol_chol, int model, int interaction, int ngpreds,
 		  bool iscox, int nullmodel=0)
     {
-	//		cout << "model = " << model << "\n";
-	//		cdata.X.print();
 	coxph_data cdata = cdatain.get_unmasked_data();
-	mematrix<double> X = t_apply_model(cdata.X,model, interaction, ngpreds, iscox, nullmodel);
+	mematrix<double> X = t_apply_model(cdata.X, model, interaction,
+					   ngpreds, iscox, nullmodel);
 	//		X.print();
 	int length_beta = X.nrow;
 	beta.reinit(length_beta,1);