[Genabel-commits] r965 - in pkg/ProbABEL: examples src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Oct 3 23:18:26 CEST 2012


Author: lckarssen
Date: 2012-10-03 23:18:26 +0200 (Wed, 03 Oct 2012)
New Revision: 965

Added:
   pkg/ProbABEL/examples/mmscore_InvSigma_aj.sex.age.small.dat
   pkg/ProbABEL/examples/mmscore_gen.small.mldose
   pkg/ProbABEL/examples/mmscore_pheno.small.PHE
Modified:
   pkg/ProbABEL/src/reg1.h
Log:
Fix bug #2295. The invvar matrix was not correctly subset leading to
rows and columns of zeroes in the 'unmasked' invvarmatrix (when using
--mmscore). If, for example, 7 IDs had an NA for a SNP than the
invvarmatrix would (correctly) be 7 rows and columns smaller. However,
additionally the final seven rows and (all) the last columns would be
exactly 0. It turns out that the wrong invvarmatrix (from rdata, the
new object, instead of rdatain) was used when copying the rows and
columns from the original invvarmatrix. 

The copyright of this fix lies with the ErasmusMC, Rotterdam.

Thanks to Maarten Kooyman for spotting this! Could this also be
related to my own finding in bug #1889 (also related to the invvar
matrix)? To be checked...


The small example files also included in this commit were used for
easy checking of this (the regular mmscore_InvSigma_aj.sex.age.dat is
500x500 elements, this one is 10x10). 


Added: pkg/ProbABEL/examples/mmscore_InvSigma_aj.sex.age.small.dat
===================================================================
--- pkg/ProbABEL/examples/mmscore_InvSigma_aj.sex.age.small.dat	                        (rev 0)
+++ pkg/ProbABEL/examples/mmscore_InvSigma_aj.sex.age.small.dat	2012-10-03 21:18:26 UTC (rev 965)
@@ -0,0 +1,10 @@
+id4 0.0182686989332702 -1.85665437227337e-05 3.53271468398681e-05 -1.74067659984958e-05 7.43017952446252e-05 -4.50104634347183e-06 -0.000124985264981788 3.3678291247208e-05 5.91968582288439e-05 -0.000243645201364358
+id10 -1.85665437227311e-05 0.0182033225536506 -1.97096076622397e-05 -5.96162563595784e-05 -5.4476919601606e-05 -9.50226926118848e-05 8.62259783546112e-05 7.3858364200391e-05 3.10831766033596e-05 -0.000130544230147021
+id25 3.53271468398634e-05 -1.97096076622424e-05 0.0181688614380818 8.52863460066612e-05 -2.43025930791677e-05 2.81566776967037e-05 4.85600196715389e-06 -6.56752842697857e-06 -1.78033754378254e-05 -1.72662880972169e-05
+id33 -1.74067659984974e-05 -5.96162563595756e-05 8.52863460066646e-05 0.0181194433068409 -5.18411780980829e-05 1.50204770703661e-05 8.56024933212197e-06 4.86573389645691e-05 3.91817984047859e-05 2.79910558167308e-05
+id35 7.43017952446215e-05 -5.44769196016117e-05 -2.43025930791681e-05 -5.18411780980874e-05 0.018106435724767 2.19056691758357e-05 -0.000128600049471664 7.17399863803957e-06 -4.24308810003088e-05 -9.37329781039098e-06
+id58 -4.50104634347325e-06 -9.50226926118893e-05 2.81566776967108e-05 1.50204770703701e-05 2.19056691758332e-05 0.0179881874374435 -1.17277182399346e-05 -4.99197036924978e-05 5.71784745599148e-05 7.8067768163972e-05
+id68 -0.000124985264981788 8.62259783546178e-05 4.85600196715395e-06 8.56024933212903e-06 -0.000128600049471660 -1.17277182399364e-05 0.0181586089581595 2.04445273957100e-05 9.39041676761655e-05 2.69629411670812e-05
+id74 3.36782912472085e-05 7.38583642003905e-05 -6.56752842698657e-06 4.86573389645756e-05 7.17399863804557e-06 -4.99197036925027e-05 2.04445273957019e-05 0.0179962093643653 -0.000143202021910123 -0.000129318612519593
+id107 5.91968582288441e-05 3.10831766033532e-05 -1.78033754378163e-05 3.91817984047789e-05 -4.2430881000305e-05 5.71784745599093e-05 9.3904167676175e-05 -0.000143202021910134 0.0182122458019134 -8.89499732951819e-05
+id119 -0.000243645201364356 -0.000130544230147016 -1.72662880972305e-05 2.79910558167332e-05 -9.37329781039559e-06 7.80677681639793e-05 2.69629411670782e-05 -0.000129318612519587 -8.89499732951819e-05 0.0183196994000799

Added: pkg/ProbABEL/examples/mmscore_gen.small.mldose
===================================================================
--- pkg/ProbABEL/examples/mmscore_gen.small.mldose	                        (rev 0)
+++ pkg/ProbABEL/examples/mmscore_gen.small.mldose	2012-10-03 21:18:26 UTC (rev 965)
@@ -0,0 +1,10 @@
+1->id4 MLDOSE 1 0 2 0 1
+1->id10 MLDOSE 2 2 1 2 0
+1->id25 MLDOSE 1 2 0 2 0
+1->id33 MLDOSE NaN 0 2 0 0
+1->id35 MLDOSE 1 0 0 0 0
+1->id58 MLDOSE 2 1 NaN 0 0
+1->id68 MLDOSE 0 1 1 0 1
+1->id74 MLDOSE 1 0 0 0 0
+1->id107 MLDOSE 1 0 NaN 1 0
+1->id119 MLDOSE 2 1 0 0 0

Added: pkg/ProbABEL/examples/mmscore_pheno.small.PHE
===================================================================
--- pkg/ProbABEL/examples/mmscore_pheno.small.PHE	                        (rev 0)
+++ pkg/ProbABEL/examples/mmscore_pheno.small.PHE	2012-10-03 21:18:26 UTC (rev 965)
@@ -0,0 +1,11 @@
+id height_residuals
+id4 -9.8130684563732
+id10 -3.47642823404311
+id25 6.26095765925731
+id33 11.3778459732314
+id35 -4.79095430879806
+id58 0.370774680307449
+id68 -8.65341409945518
+id74 4.95199124661318
+id107 -9.8280719102797
+id119 -2.05732829257249

Modified: pkg/ProbABEL/src/reg1.h
===================================================================
--- pkg/ProbABEL/src/reg1.h	2012-09-30 14:54:31 UTC (rev 964)
+++ pkg/ProbABEL/src/reg1.h	2012-10-03 21:18:26 UTC (rev 965)
@@ -360,21 +360,25 @@
 	// model should come here
 	regdata rdata = rdatain.get_unmasked_data();
 
+	// Create a new invvarmatrix that only contains the values for
+	// individuals that are not masked.
 	mematrix<double> invvarmatrix;
 	if(invvarmatrixin.nrow!=0 && invvarmatrixin.ncol!=0)
 	{
 	    invvarmatrix.reinit(rdata.nids,rdata.nids);
 	    int i1=0, j1=0;
-	    for (int i=0;i<rdata.nids;i++)
+	    for (int i=0; i<rdatain.nids; i++) {
 		if (rdatain.masked_data[i]==0) {
 		    j1 = 0;
-		    for (int j=0;j<rdata.nids;j++) if (rdatain.masked_data[j]==0) {
+		    for (int j=0; j<rdatain.nids; j++) {
+			if (rdatain.masked_data[j]==0) {
 			    invvarmatrix.put(invvarmatrixin.get(i,j),i1,j1);
 			    j1++;
 			}
+		    }
 		    i1++;
 		}
-
+	    }
 	}
 
 	//fprintf(stdout,"estimate: %i %i %i %i\n",rdata.nids,(rdata.X).nrow,(rdata.X).ncol,(rdata.Y).ncol);



More information about the Genabel-commits mailing list