View Javadoc
1   /*
2    * The baseCode project
3    *
4    * Copyright (c) 2011 University of British Columbia
5    *
6    * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with
7    * the License. You may obtain a copy of the License at
8    *
9    * http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on
12   * an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the
13   * specific language governing permissions and limitations under the License.
14   */
15  package ubic.basecode.math.linearmodels;
16  
17  import cern.colt.list.DoubleArrayList;
18  import cern.colt.matrix.DoubleMatrix1D;
19  import cern.colt.matrix.DoubleMatrix2D;
20  import cern.colt.matrix.impl.DenseDoubleMatrix2D;
21  import cern.colt.matrix.linalg.Algebra;
22  import cern.jet.math.Functions;
23  import org.apache.commons.math3.distribution.FDistribution;
24  import org.junit.Test;
25  import ubic.basecode.dataStructure.matrix.*;
26  import ubic.basecode.io.reader.DoubleMatrixReader;
27  import ubic.basecode.io.reader.StringMatrixReader;
28  import ubic.basecode.math.DescriptiveWithMissing;
29  import ubic.basecode.math.MatrixStats;
30  
31  import java.util.List;
32  import java.util.Map;
33  import java.util.zip.GZIPInputStream;
34  
35  import static org.junit.Assert.*;
36  
37  /**
38   * @author paul
39   */
40  public class LeastSquaresFitTest {
41  
42      /**
43       * @throws Exception
44       */
45      @Test
46      public void testLSFOneContinuousWithMissing3() throws Exception {
47          DoubleMatrixReader f = new DoubleMatrixReader();
48          DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
49                  "/data/example.madata.withmissing.small.txt" ) );
50  
51          ObjectMatrix<String, String, Object> design = new ObjectMatrixImpl<>( 9, 1 );
52          design.set( 0, 0, 0.12 );
53          design.set( 1, 0, 0.24 );
54          design.set( 2, 0, 0.48 );
55          design.set( 3, 0, 0.96 );
56          design.set( 4, 0, 0.12 );
57          design.set( 5, 0, 0.24 );
58          design.set( 6, 0, 0.48 );
59          design.set( 7, 0, 0.96 );
60          design.set( 8, 0, 0.96 );
61  
62          design.setRowNames( testMatrix.getColNames() );
63          design.addColumnName( "Value" );
64  
65          LeastSquaresFit fit = new LeastSquaresFit( design, testMatrix );
66  
67          Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
68  
69          LinearModelSummary s = sums.get( "228980_at" ); // has missing
70          assertNotNull( s );
71          assertEquals( 0.1495, s.getF(), 0.01 );
72          assertEquals( 0.7123, s.getP(), 0.001 );
73          assertEquals( 10.9180, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
74          assertEquals( 0.712, s.getContrastCoefficients().get( 1, 3 ), 0.001 ); // pvalue
75          assertEquals( 6, s.getResidualDof().intValue() );
76          GenericAnovaResult a = s.getAnova();
77          assertEquals( 0.1495, a.getMainEffectF( "Value" ), 0.0001 );
78          assertEquals( 1, a.getMainEffectDof( "Value" ).intValue() );
79          assertEquals( 6, a.getResidualDf().intValue() );
80  
81          FDistribution fd = new FDistribution( 1, 6 );
82          double p = 1.0 - fd.cumulativeProbability( 0.1495 );
83          assertEquals( 0.7123, p, 0.0001 );
84          assertEquals( 0.7123, a.getMainEffectP( "Value" ), 0.0001 );
85  
86          s = sums.get( "1553129_at" );
87          assertNotNull( s );
88          assertEquals( 2.095, s.getF(), 0.01 );
89          assertEquals( 0.1911, s.getP(), 0.001 );
90          assertEquals( 3.78719, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
91          assertEquals( 0.191, s.getContrastCoefficients().get( 1, 3 ), 0.001 ); // this ordering might change?
92          a = s.getAnova();
93          assertNotNull( a );
94          assertEquals( 0.1911, a.getMainEffectP( "Value" ), 0.0001 );
95  
96          s = fit.summarize( 14 );
97          assertNotNull( s );
98          assertEquals( "214502_at", s.getKey() );// has missing
99          assertEquals( 1.992, s.getF(), 0.01 );
100         assertEquals( 0.2172, s.getP(), 0.001 );
101         assertEquals( 4.2871, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
102         assertEquals( 0.217, s.getContrastCoefficients().get( 1, 3 ), 0.001 );
103 
104         s = sums.get( "232018_at" );
105         assertNotNull( s );
106         assertEquals( 1.381, s.getF(), 0.01 );
107         assertEquals( 0.2783, s.getP(), 0.001 );
108         assertEquals( 6.6537, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
109         assertEquals( 0.278, s.getContrastCoefficients().get( 1, 3 ), 0.001 );
110         a = s.getAnova();
111         assertNotNull( a );
112         assertEquals( 0.2783, a.getMainEffectP( "Value" ), 0.0001 );
113         s = sums.get( "228980_at" ); // has missing
114         assertNotNull( s );
115         assertEquals( 0.1495, s.getF(), 0.01 );
116         assertEquals( 0.7123, s.getP(), 0.001 );
117         assertEquals( 10.9180, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
118         assertEquals( 0.712, s.getContrastCoefficients().get( 1, 3 ), 0.001 );
119         a = s.getAnova();
120         assertNotNull( a );
121         assertEquals( 0.7123, a.getMainEffectP( "Value" ), 0.0001 );
122     }
123 
124     /**
125      * @throws Exception
126      */
127     @Test
128     public void testLSFThreeLevelsOnecontinous2() throws Exception {
129         DoubleMatrixReader f = new DoubleMatrixReader();
130         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
131                 "/data/example.madata.small.txt" ) );
132 
133         ObjectMatrix<String, String, Object> design = new ObjectMatrixImpl<>( 9, 2 );
134 
135         design.set( 0, 0, "A" );
136         design.set( 1, 0, "A" );
137         design.set( 2, 0, "A" );
138         design.set( 3, 0, "B" );
139         design.set( 4, 0, "B" );
140         design.set( 5, 0, "B" );
141         design.set( 6, 0, "C" );
142         design.set( 7, 0, "C" );
143         design.set( 8, 0, "C" );
144         design.set( 0, 1, 0.12 );
145         design.set( 1, 1, 0.24 );
146         design.set( 2, 1, 0.48 );
147         design.set( 3, 1, 0.96 );
148         design.set( 4, 1, 0.12 );
149         design.set( 5, 1, 0.24 );
150         design.set( 6, 1, 0.48 );
151         design.set( 7, 1, 0.96 );
152         design.set( 8, 1, 0.96 );
153 
154         design.setRowNames( testMatrix.getColNames() );
155         design.addColumnName( "Factor" );
156         design.addColumnName( "Value" );
157 
158         LeastSquaresFit fit = new LeastSquaresFit( design, testMatrix );
159 
160         DoubleMatrix2D coeffs = fit.getCoefficients();
161         assertEquals( 3.7458080, coeffs.get( 0, 0 ), 0.0001 );
162         assertEquals( -0.4388889, coeffs.get( 1, 2 ), 0.0001 );
163         assertEquals( 0.5709091, coeffs.get( 2, 10 ), 0.0001 );
164         assertEquals( 0.04856061, coeffs.get( 2, 18 ), 0.0001 );
165         assertEquals( -1.1363636, coeffs.get( 3, 10 ), 0.0001 );
166         assertEquals( 0.11174242, coeffs.get( 3, 18 ), 0.0001 );
167 
168         DoubleMatrix2D fitted = fit.getFitted();
169 
170         assertEquals( 3.764747, fitted.get( 0, 0 ), 0.0001 );
171         assertEquals( 6.043990, fitted.get( 1, 3 ), 0.0001 );
172         assertEquals( 10.858586, fitted.get( 7, 2 ), 0.0001 );
173         assertEquals( 6.307879, fitted.get( 18, 8 ), 0.0001 );
174 
175         List<GenericAnovaResult> anova = fit.anova();
176         assertEquals( 19, anova.size() );
177     }
178 
179     /**
180      * @throws Exception
181      */
182     @Test
183     public void testLSFThreeLevelsOneContinuousWithMissing3() throws Exception {
184         DoubleMatrixReader f = new DoubleMatrixReader();
185         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
186                 "/data/example.madata.withmissing.small.txt" ) );
187 
188         ObjectMatrix<String, String, Object> design = new ObjectMatrixImpl<>( 9, 2 );
189 
190         design.set( 0, 0, "A" );
191         design.set( 1, 0, "A" );
192         design.set( 2, 0, "A" );
193         design.set( 3, 0, "B" );
194         design.set( 4, 0, "B" );
195         design.set( 5, 0, "B" );
196         design.set( 6, 0, "C" );
197         design.set( 7, 0, "C" );
198         design.set( 8, 0, "C" );
199         design.set( 0, 1, 0.12 );
200         design.set( 1, 1, 0.24 );
201         design.set( 2, 1, 0.48 );
202         design.set( 3, 1, 0.96 );
203         design.set( 4, 1, 0.12 );
204         design.set( 5, 1, 0.24 );
205         design.set( 6, 1, 0.48 );
206         design.set( 7, 1, 0.96 );
207         design.set( 8, 1, 0.96 );
208 
209         design.setRowNames( testMatrix.getColNames() );
210         design.addColumnName( "Factor" );
211         design.addColumnName( "Value" );
212 
213         LeastSquaresFit fit = new LeastSquaresFit( design, new DenseDoubleMatrix2D( testMatrix.asArray() ) );
214 
215         DoubleMatrix2D coeffs = fit.getCoefficients();
216         assertEquals( 3.7458080, coeffs.get( 0, 0 ), 0.0001 );
217         assertEquals( -0.4388889, coeffs.get( 1, 2 ), 0.0001 );
218         assertEquals( 0.5709091, coeffs.get( 2, 10 ), 0.0001 );
219         assertEquals( 0.04856061, coeffs.get( 2, 18 ), 0.0001 );
220         assertEquals( -1.1363636, coeffs.get( 3, 10 ), 0.0001 );
221         assertEquals( 0.11174242, coeffs.get( 3, 18 ), 0.0001 );
222 
223         DoubleMatrix2D fitted = fit.getFitted();
224 
225         assertEquals( 3.764747, fitted.get( 0, 0 ), 0.0001 );
226         assertEquals( 6.043990, fitted.get( 1, 3 ), 0.0001 );
227         assertEquals( 10.8333, fitted.get( 7, 2 ), 0.0001 );
228         assertEquals( 6.307879, fitted.get( 18, 8 ), 0.0001 );
229 
230         List<GenericAnovaResult> anova = fit.anova();
231 
232         assertEquals( 19, anova.size() );
233     }
234 
235     /**
236      * @throws Exception
237      */
238     @Test
239     public void testLSFTwoLevels() throws Exception {
240         DoubleMatrixReader f = new DoubleMatrixReader();
241         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
242                 "/data/example.madata.small.txt" ) );
243 
244         ObjectMatrix<String, String, Object> design = new ObjectMatrixImpl<>( 9, 1 );
245 
246         design.set( 0, 0, "A" );
247         design.set( 1, 0, "A" );
248         design.set( 2, 0, "A" );
249         design.set( 3, 0, "A" );
250         design.set( 4, 0, "B" );
251         design.set( 5, 0, "B" );
252         design.set( 6, 0, "B" );
253         design.set( 7, 0, "B" );
254         design.set( 8, 0, "B" );
255         design.setRowNames( testMatrix.getColNames() );
256         design.addColumnName( "Factor" );
257 
258         LeastSquaresFit fit = new LeastSquaresFit( design, testMatrix );
259 
260         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
261 
262         LinearModelSummary s = sums.get( "1553129_at" );
263         assertEquals( 0.182999, s.getF(), 0.0001 );
264         assertEquals( 0.6817, s.getP(), 0.001 );
265         assertEquals( 3.84250, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
266         assertEquals( 0.682, s.getContrastCoefficients().get( 1, 3 ), 0.001 ); // pvalue.
267 
268         Double[] effects = s.getEffects();
269         assertEquals( -11.58333, effects[0], 0.0001 ); // hm.
270         assertEquals( 0.04999, effects[1], 0.0001 ); //hm
271 
272         Double[] stdevUnscaled = s.getStdevUnscaled(); //
273         assertEquals( 0.5, stdevUnscaled[0], 0.0001 );
274         assertEquals( 0.6708203932, stdevUnscaled[1], 0.0001 );
275 
276         Double sigma = s.getSigma();
277         assertEquals( 0.11673841331, sigma, 0.0001 );
278 
279         s = sums.get( "232018_at" );
280         assertEquals( -18.9866667, s.getEffects()[0], 0.0001 ); // hm.
281         assertEquals( 0.1714319, s.getEffects()[1], 0.0001 ); //hm
282         assertEquals( 0.07879, s.getF(), 0.01 );
283         assertEquals( 0.787, s.getP(), 0.001 );
284         assertEquals( 6.2650, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
285         assertEquals( 0.787, s.getContrastCoefficients().get( 1, 3 ), 0.001 );
286         sigma = s.getSigma();
287         assertEquals( 0.61072556381, sigma, 0.0001 );
288     }
289 
290     /**
291      * Many missing values; Two factors, two levels + interaction.
292      *
293      * @throws Exception
294      */
295     @Test
296     public void testLSFTwoLevels2() throws Exception {
297         DoubleMatrixReader f = new DoubleMatrixReader();
298         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
299                 "/data/lmtest1.dat.manymissing.txt" ) );
300 
301         StringMatrixReader of = new StringMatrixReader();
302         StringMatrix<String, String> sampleInfo = of.read( this.getClass()
303                 .getResourceAsStream( "/data/lmtest1.des.txt" ) );
304 
305         DesignMatrix d = new DesignMatrix( sampleInfo, true );
306         d.addInteraction();
307 
308         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix );
309         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
310         assertEquals( 100, sums.size() );
311 
312         for ( LinearModelSummary lms : sums.values() ) {
313             GenericAnovaResult a = lms.getAnova();
314             assertNotNull( a );
315             Double interactionEffectP = a.getInteractionEffectP();
316             assertNotNull( interactionEffectP );
317         }
318 
319     }
320 
321     /**
322      * @throws Exception
323      */
324     @Test
325     public void testLSFTwoLevels3() throws Exception {
326         DoubleMatrixReader f = new DoubleMatrixReader();
327         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
328                 "/data/anova-test-data.txt" ) );
329 
330         StringMatrixReader of = new StringMatrixReader();
331         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
332                 "/data/anova-test-des.txt" ) );
333 
334         DesignMatrix d = new DesignMatrix( sampleInfo, true );
335         d.addInteraction( "factor1", "factor2" );
336 
337         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix );
338         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
339         assertEquals( 100, sums.size() );
340 
341         for ( LinearModelSummary lms : sums.values() ) {
342             GenericAnovaResult a = lms.getAnova();
343             assertNotNull( a );
344             Double interactionEffectP = a.getInteractionEffectP();
345             assertNotNull( interactionEffectP );
346         }
347 
348         assertEquals( 0.0048, sums.get( "probe_4" ).getMainEffectP( "factor1" ), 0.0001 );
349         assertEquals( 5.158e-10, sums.get( "probe_10" ).getMainEffectP( "factor1" ), 1e-12 );
350         assertEquals( 0.6888, sums.get( "probe_98" ).getMainEffectP( "factor2" ), 1e-4 );
351         assertEquals( 0.07970, sums.get( "probe_10" ).getMainEffectP( "factor2" ), 1e-4 );
352 
353     }
354 
355     /**
356      * @throws Exception
357      */
358     @Test
359     public void testLSFTwoLevelsOneContinuous() throws Exception {
360         DoubleMatrixReader f = new DoubleMatrixReader();
361         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
362                 "/data/example.madata.small.txt" ) );
363 
364         ObjectMatrix<String, String, Object> design = new ObjectMatrixImpl<>( 9, 2 );
365 
366         design.set( 0, 0, "A" );
367         design.set( 1, 0, "A" );
368         design.set( 2, 0, "A" );
369         design.set( 3, 0, "A" );
370         design.set( 4, 0, "B" );
371         design.set( 5, 0, "B" );
372         design.set( 6, 0, "B" );
373         design.set( 7, 0, "B" );
374         design.set( 8, 0, "B" );
375         design.set( 0, 1, 0.12 );
376         design.set( 1, 1, 0.24 );
377         design.set( 2, 1, 0.48 );
378         design.set( 3, 1, 0.96 );
379         design.set( 4, 1, 0.12 );
380         design.set( 5, 1, 0.24 );
381         design.set( 6, 1, 0.48 );
382         design.set( 7, 1, 0.96 );
383         design.set( 8, 1, 0.96 );
384 
385         design.setRowNames( testMatrix.getColNames() );
386         design.addColumnName( "Factor" );
387         design.addColumnName( "Value" );
388 
389         LeastSquaresFit fit = new LeastSquaresFit( design, testMatrix );
390 
391         DoubleMatrix2D coeffs = fit.getCoefficients();
392         assertEquals( 3.77868, coeffs.get( 0, 0 ), 0.0001 );
393         assertEquals( 0.24476, coeffs.get( 1, 2 ), 0.0001 );
394         assertEquals( -0.680449, coeffs.get( 2, 10 ), 0.0001 );
395         assertEquals( 0.114084, coeffs.get( 2, 18 ), 0.0001 );
396 
397         DoubleMatrix2D fitted = fit.getFitted();
398 
399         assertEquals( 3.795698, fitted.get( 0, 0 ), 0.0001 );
400         assertEquals( 5.497165, fitted.get( 1, 3 ), 0.0001 );
401         assertEquals( 10.879917, fitted.get( 7, 2 ), 0.0001 );
402         assertEquals( 6.346546, fitted.get( 18, 8 ), 0.0001 );
403 
404         List<GenericAnovaResult> anova = fit.anova();
405         assertEquals( 19, anova.size() );
406         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
407 
408         LinearModelSummary s = sums.get( "1553129_at" );
409         assertEquals( 0.9389, s.getF(), 0.01 );
410         assertEquals( 0.4418, s.getP(), 0.001 );
411         assertEquals( 3.77868, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
412         assertEquals( 0.810, s.getContrastCoefficients().get( 1, 3 ), 0.001 ); // this ordering might change?
413         GenericAnovaResult a = s.getAnova();
414         assertEquals( 0.2429, a.getMainEffectP( "Value" ), 0.0001 );
415 
416         s = sums.get( "232018_at" );
417         assertEquals( 0.7167, s.getF(), 0.01 );
418         assertEquals( 0.5259, s.getP(), 0.001 );
419         assertEquals( 6.5712, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
420         assertEquals( 0.664, s.getContrastCoefficients().get( 1, 3 ), 0.001 );
421         a = s.getAnova();
422         assertEquals( 0.2893, a.getMainEffectP( "Value" ), 0.0001 );
423 
424         // 232018_at // based on rstudent() in R.
425         DoubleMatrix2D studentizedResiduals = fit.getStudentizedResiduals();
426         // log.info( studentizedResiduals.viewRow( 10 ) );
427         double[] expectedStudentizedResiduals = new double[] { -0.34655041, 1.46251738, -0.61403124, -0.34663812,
428                 -1.51245468, 0.06875469, 1.45818880, 1.02811044, -1.31696150 };
429 
430         for ( int i = 0; i < 9; i++ ) {
431             assertEquals( expectedStudentizedResiduals[i], studentizedResiduals.viewRow( 10 ).get( i ), 0.001 );
432         }
433 
434         // assertEquals( 1.1, DescriptiveWithMissing.variance( new DoubleArrayList( expectedStudentizedResiduals ) ),
435         // 0.1 );
436 
437         assertEquals( 1.1,
438                 DescriptiveWithMissing.variance( new DoubleArrayList( studentizedResiduals.viewRow( 10 ).toArray() ) ),
439                 0.1 );
440 
441         // 1553129_at // based on rstudent() in R.
442         studentizedResiduals = fit.getStudentizedResiduals();
443         expectedStudentizedResiduals = new double[] { 0.46128657, -5.49429390, 0.84157385, 1.10053286, 1.10538546,
444                 -0.01706794, -0.05318259, -0.56926585, -0.35107932 };
445         for ( int i = 0; i < 9; i++ ) {
446             assertEquals( expectedStudentizedResiduals[i], studentizedResiduals.viewRow( 0 ).get( i ), 0.001 );
447         }
448 
449     }
450 
451     /**
452      * Weighted least squares test for 2D matrices
453      *
454      * @throws Exception
455      */
456     @Test
457     public void testMatrixWeightedRegress() throws Exception {
458         DoubleMatrix2D dat = new DenseDoubleMatrix2D( new double[][] { { 1, 2, 3, 4, 5 }, { 1, 1, 6, 3, 2 } } );
459         DoubleMatrix2D des = new DenseDoubleMatrix2D( new double[][] { { 1, 1, 1, 1, 1 }, { 1, 2, 2, 3, 3 },
460                 { 2, 1, 5, 3, 4 } } );
461 
462         DoubleMatrix2D w = dat.copy();
463         w.assign( dat );
464         Algebra solver = new Algebra();
465         des = solver.transpose( des );
466         LeastSquaresFit fit = new LeastSquaresFit( des, dat, w );
467 
468         /*
469          * TODO R code please
470          */
471 
472         // FIXME why is precision of these tests so low? It was 0.1! I changed it to 0.001
473 
474         // coefficients
475         DoubleMatrix2D actuals = solver.transpose( fit.getCoefficients() );
476         double[][] expected = new double[][] { { -1.7070, 1.7110, 0.3054 }, { 0.2092, -0.6642, 1.3640 } };
477         for ( int i = 0; i < expected.length; i++ ) {
478             assertArrayEquals( expected[i], actuals.viewRow( i ).toArray(), 0.001 );
479         }
480 
481         // fitted
482         actuals = fit.getFitted();
483         expected = new double[][] { { 0.6151, 2.0210, 3.2430, 4.3430, 4.6490 }, { 2.273, 0.245, 5.701, 2.309, 3.673 } };
484         for ( int i = 0; i < expected.length; i++ ) {
485             assertArrayEquals( expected[i], actuals.viewRow( i ).toArray(), 0.001 );
486         }
487 
488         // residuals
489         actuals = fit.getResiduals();
490         expected = new double[][] { { 0.38490, -0.02092, -0.24270, -0.34310, 0.35150 },
491                 { -1.2730, 0.7550, 0.2986, 0.6910, -1.6730 } };
492         for ( int i = 0; i < expected.length; i++ ) {
493             assertArrayEquals( expected[i], actuals.viewRow( i ).toArray(), 0.001 );
494         }
495 
496     }
497 
498     /**
499      * Has a lot of missing values.
500      *
501      * @throws Exception
502      */
503     @Test
504     public void testOneWayAnova() throws Exception {
505         DoubleMatrixReader f = new DoubleMatrixReader();
506         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
507                 "/data/anova-test-data.txt" ) );
508 
509         ObjectMatrix<String, String, Object> design = new ObjectMatrixImpl<>( 8, 1 );
510         for ( int i = 0; i < 8; i++ ) {
511             design.set( i, 0, "A" + i % 3 );
512         }
513         design.addColumnName( "Factor1" );
514 
515         DesignMatrix d = new DesignMatrix( design, true );
516 
517         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix );
518         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
519         assertEquals( 100, sums.size() );
520 
521         for ( LinearModelSummary lms : sums.values() ) {
522             GenericAnovaResult a = lms.getAnova();
523             assertNotNull( a );
524             Double interactionEffectP = a.getInteractionEffectP();
525             assertNull( interactionEffectP );
526         }
527 
528         LinearModelSummary sum4 = sums.get( "probe_4" );
529         assertNotNull( sum4.getContrastCoefficients() );
530         assertEquals( 0.6531, sum4.getP(), 0.0001 );
531         assertEquals( 0.2735, sum4.getF(), 0.0001 );
532         assertEquals( 0.2735, sum4.getAnova().getMainEffectF( "Factor1" ), 0.0001 );
533         assertEquals( 2, sum4.getAnova().getResidualDf().intValue() );
534         assertEquals( 1, sum4.getAnova().getMainEffectDof( "Factor1" ).intValue() );
535         assertEquals( 0.6531, sum4.getMainEffectP( "Factor1" ), 0.0001 );
536 
537         LinearModelSummary sum21 = sums.get( "probe_21" );
538         assertNotNull( sum21.getContrastCoefficients() );
539         assertEquals( 0.6492, sum21.getP(), 0.0001 );
540         assertEquals( 0.4821, sum21.getF(), 0.0001 );
541         assertEquals( 0.4821, sum21.getAnova().getMainEffectF( "Factor1" ), 0.0001 );
542         assertEquals( 4, sum21.getAnova().getResidualDf().intValue() );
543         assertEquals( 2, sum21.getAnova().getMainEffectDof( "Factor1" ).intValue() );
544         assertEquals( 0.6492, sum21.getMainEffectP( "Factor1" ), 0.0001 );
545 
546         LinearModelSummary sum98 = sums.get( "probe_98" );
547         assertNotNull( sum98.getContrastCoefficients() );
548         assertEquals( 0.1604, sum98.getP(), 0.0001 );
549         assertEquals( 2.993, sum98.getF(), 0.0001 );
550         assertEquals( 4, sum98.getAnova().getResidualDf().intValue() );
551         assertEquals( 2, sum98.getAnova().getMainEffectDof( "Factor1" ).intValue() );
552         assertEquals( 2.9931, sum98.getAnova().getMainEffectF( "Factor1" ).doubleValue(), 0.0001 );
553         assertEquals( 0.1604, sum98.getMainEffectP( "Factor1" ), 1e-4 );
554 
555         LinearModelSummary sum10 = sums.get( "probe_10" );
556         assertNotNull( sum10.getContrastCoefficients() );
557         assertEquals( 0.8014, sum10.getP(), 0.0001 );
558         assertEquals( 0.2314, sum10.getF(), 0.0001 );
559         assertEquals( 5, sum10.getAnova().getResidualDf().intValue() );
560         assertEquals( 2, sum10.getAnova().getMainEffectDof( "Factor1" ).intValue() );
561         assertEquals( 0.8014, sum10.getMainEffectP( "Factor1" ), 1e-4 );
562 
563         /*
564          * Painful case follows. Not full rank due to missing values:
565          */
566         LinearModelSummary sum60 = sums.get( "probe_60" );
567         assertNotNull( sum60.getContrastCoefficients() );
568 
569         /*
570          * See lmtests.R
571          * 
572          * Note that using lmFit and lm give different results because (apparently) of the difference in strategy
573          * dealing
574          * with missing values. In lm, the factors are subset; in lmFit, the design matrix is subset. We do the latter.
575          *
576          * lmFit(dat, model.matrix(~ factor3))
577          * 
578          * vs
579          * 
580          * lm(t(dat["probe_60",]) ~ factor3)
581          */
582         assertEquals( 9.00145, sum60.getContrastCoefficients().get( 0, 0 ), 0.0001 ); // same as lmFit; R lm gives 8.9913; lm.fit gives tehse values.
583         assertEquals( -0.01020, sum60.getContrastCoefficients().get( 1, 0 ), 0.0001 ); // same as lmFit; R lm gives +0.0102
584         assertEquals( Double.NaN, sum60.getContrastCoefficients().get( 2, 0 ), 0.0001 ); // good ...
585 
586         /*
587          * This is the value of R-squared from lm(). lmFit doesn't give us this, but...
588          */
589         //   assertEquals( -0.4996, sum60.getAdjRSquared(), 1e-5 );
590 
591         // using results of lm()
592         // anova(object)
593         //        Analysis of Variance Table
594         //
595         //        Response: t(dat["probe_60", ])
596         //                  Df  Sum Sq  Mean Sq F value Pr(>F)
597         //        factor3    1 0.00010 0.000104   5e-04 0.9846 <- not quite
598         //        Residuals  2 0.44135 0.220675   <- we get this
599 
600         assertEquals( 2, sum60.getResidualDof().intValue() );
601         assertEquals( 1, sum60.getNumeratorDof().intValue() );
602 
603         assertEquals( 2, sum60.getAnova().getResidualDf().intValue() );
604         assertEquals( 1, sum60.getAnova().getMainEffectDof( "Factor1" ).intValue() );
605 
606         // The value lm gives is 0.004715; for limma, topTableF gives NA and no pvalue. 
607         // However, for this case, this is actually the t-statistic (two groups) and the right value is 0.022 (from lm)
608         assertEquals( 0.0004715, sum60.getF(), 1e-7 ); // on 1 and 2 dof,  p-value: 0.9846
609         assertEquals( 0.9846482, sum60.getP(), 1e-5 );
610 
611         //
612         assertEquals( 0.9846482, sum60.getMainEffectP( "Factor1" ), 1e-5 );
613 
614         // Result for topTable(eBayes(fit), coef=2)
615         //                logFC    AveExpr            t      P.Value   adj.P.Val          B
616         // probe_60  -0.01020000   8.996350  -0.02338077 9.831252e-01 0.983125229 -7.4711642
617 
618         // topTableF(eBayes(o2), number = 100)
619         //          X.Intercept.     factor3v factor3w_base    AveExpr            F      P.Value    adj.P.Val
620         // probe_60     9.001450  -0.01020000            NA   8.996350           NA           NA           NA
621 
622         // Using lm(formula = t(dat["probe_60", ]) ~ factor3)
623 
624         //        Call:
625         //            lm(formula = t(dat["probe_60", ]) ~ factor3)
626         //
627         //            Residuals:
628         //            X0b.bioassay.0b. X1a.bioassay.1a. X2a.bioassay.2a. X2b.bioassay.2b. 
629         //                     -0.1209          -0.4539           0.1209           0.4539 
630         //
631         //            Coefficients:
632         //                          Estimate Std. Error t value Pr(>|t|)   
633         //            (Intercept)     8.9913     0.3322  27.068  0.00136 **
634         //            factor3w_base   0.0102     0.4698   0.022  0.98465   
635         //            ---
636         //            Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
637         //
638         //            Residual standard error: 0.4698 on 2 degrees of freedom
639         //              (4 observations deleted due to missingness)
640         //            Multiple R-squared:  0.0002357, Adjusted R-squared:  -0.4996 
641         //            F-statistic: 0.0004715 on 1 and 2 DF,  p-value: 0.9846
642 
643         // results of using      lm.fit(X,y)
644         /// > X
645         //        (Intercept) factor3v factor3w_base
646         //        2           1        1             0
647         //        3           1        0             1
648         //        5           1        1             0
649         //        6           1        0             1
650         //        > y
651         //        [1] 8.8704 8.5475 9.1121 9.4554
652         //        > lm.fit(X,y)
653         //        $coefficients
654         //          (Intercept)      factor3v factor3w_base 
655         //              9.00145      -0.01020            NA 
656         //
657         //        $residuals
658         //        [1] -0.12085 -0.45395  0.12085  0.45395
659         //
660         //        $effects
661         //        (Intercept)    factor3v                         
662         //           -17.9927     -0.0102      0.0784      0.6597 
663         //
664         //        $rank
665         //        [1] 2
666         //
667         //        $fitted.values
668         //        [1] 8.99125 9.00145 8.99125 9.00145
669         //
670         //        $assign
671         //        NULL
672         //
673         //        $qr
674         //        $qr
675         //          (Intercept)   factor3v factor3w_base
676         //        2        -2.0 -1.0000000 -1.000000e+00
677         //        3         0.5  1.0000000 -1.000000e+00
678         //        5         0.5 -0.3333333 -1.110223e-16
679         //        6         0.5  0.6666667  0.000000e+00
680         //
681         //        $qraux
682         //        [1] 1.500000 1.666667 2.000000
683         //
684         //        $pivot
685         //        [1] 1 2 3
686         //
687         //        $tol
688         //        [1] 1e-07
689         //
690         //        $rank
691         //        [1] 2
692         //
693         //        attr(,"class")
694         //        [1] "qr"
695         //
696         //        $df.residual
697         //        [1] 2
698 
699     }
700 
701     /**
702      * @throws Exception
703      */
704     @Test
705     public void testSingular() throws Exception {
706         DoubleMatrixReader f = new DoubleMatrixReader();
707         DoubleMatrix<String, String> testMatrix = f
708                 .read( this.getClass().getResourceAsStream( "/data/lmtest2.dat.txt" ) );
709 
710         StringMatrixReader of = new StringMatrixReader();
711         StringMatrix<String, String> sampleInfo = of.read( this.getClass()
712                 .getResourceAsStream( "/data/lmtest2.des.txt" ) );
713         DesignMatrix d = new DesignMatrix( sampleInfo, true );
714 
715         assertEquals( 9, d.getMatrix().columns() );
716 
717         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix );
718 
719         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
720         assertEquals( 81, sums.size() );
721 
722         for ( LinearModelSummary lms : sums.values() ) {
723             GenericAnovaResult a = lms.getAnova();
724             assertNotNull( a );
725         }
726 
727         LinearModelSummary s = sums.get( "A01157cds_s_at" );
728 
729         assertNotNull( s.getContrastCoefficients() );
730 
731         assertEquals( 7.3740000, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
732         // assertEquals( 0.1147667, s.getCoefficients().get( 1, 3 ), 0.001 );
733         assertEquals( 6, s.getResidualDof().intValue() );
734         assertEquals( 7, s.getNumeratorDof().intValue() );
735         assertEquals( 0.8634, s.getF(), 0.01 );
736         assertEquals( 0.5795, s.getP(), 0.001 );
737 
738     }
739 
740     /**
741      * Originally causes failures during summarization step. There are two pivoted columns.
742      *
743      * @throws Exception
744      */
745     @Test
746     public void testSingular2() throws Exception {
747         DoubleMatrixReader f = new DoubleMatrixReader();
748         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
749                 "/data/1027_GSE6189.data.test.txt" ) );
750 
751         StringMatrixReader of = new StringMatrixReader();
752         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
753                 "/data/1027_GSE6189_expdesign.data.txt" ) );
754         DesignMatrix d = new DesignMatrix( sampleInfo, true );
755 
756         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix.getRowRange( 0, 0 ) );
757 
758         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
759         assertEquals( 1, sums.size() );
760         for ( LinearModelSummary lms : sums.values() ) {
761             GenericAnovaResult a = lms.getAnova();
762             assertNotNull( a );
763         }
764         LinearModelSummary s = sums.get( "1367452_at" );
765         assertNotNull( s );
766         assertNotNull( s.getContrastCoefficients() );
767 
768         // log.info( s.getCoefficients() );
769 
770         // log.info( s.getAnova() );
771 
772         // our model matrix ends up with different coefficients than R which are
773         // double[] rcoef = new double[] { 14.96355, 0.14421, -0.11525, 0.24257, Double.NaN, 0.04093, 0.06660,
774         // Double.NaN };
775 
776         // here are the coefs we get in R if we use the exact model matrix we get drom our DesignMatrix
777         double[] coef = new double[] { 15.10776244, -0.01689300, 0.09835841, -0.20163964, Double.NaN, -0.04092962,
778                 Double.NaN, 0.06660370 };
779 
780         for ( int i = 0; i < s.getContrastCoefficients().rows(); i++ ) {
781             assertEquals( coef[i], s.getContrastCoefficients().get( i, 0 ), 0.0001 );
782         }
783 
784     }
785 
786     @Test
787     public void testThreeWaySingular() throws Exception {
788         DoubleMatrixReader f = new DoubleMatrixReader();
789         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
790                 "/data/1064_GSE7863.data.test.txt" ) );
791 
792         StringMatrixReader of = new StringMatrixReader();
793         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
794                 "/data/1064_GSE7863_expdesign.data.test.txt" ) );
795         DesignMatrix d = new DesignMatrix( sampleInfo, true );
796 
797         assertEquals( 5, d.getMatrix().columns() );
798 
799         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix );
800 
801         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
802         assertEquals( 416, sums.size() );
803 
804         for ( LinearModelSummary lms : sums.values() ) {
805             GenericAnovaResult a = lms.getAnova();
806             assertNotNull( a );
807         }
808 
809         LinearModelSummary s = sums.get( "1415696_at" );
810 
811         assertNotNull( s.getContrastCoefficients() );
812         assertEquals( 0.000794, s.getContrastCoefficients().get( 2, 3 ), 0.001 );
813         assertEquals( 11, s.getResidualDof().intValue() );
814         assertEquals( 4, s.getNumeratorDof().intValue() );
815         assertEquals( 24.38, s.getF(), 0.01 );
816         assertEquals( 2.025e-05, s.getP(), 0.001 );
817         GenericAnovaResult anova = s.getAnova();
818         assertEquals( 29.0386, anova.getMainEffectF( "Treatment" ), 0.0001 );
819 
820         s = sums.get( "1415837_at" );
821         assertNotNull( s.getContrastCoefficients() );
822         assertEquals( 11, s.getResidualDof().intValue() );
823         assertEquals( 4, s.getNumeratorDof().intValue() );
824         assertEquals( 22.72, s.getF(), 0.01 );
825         assertEquals( 2.847e-05, s.getP(), 0.001 );
826         anova = s.getAnova();
827         assertEquals( 6.5977, anova.getMainEffectF( "Treatment" ), 0.0001 );
828 
829         s = sums.get( "1416179_a_at" );
830         assertNotNull( s.getContrastCoefficients() );
831         assertEquals( 11, s.getResidualDof().intValue() );
832         assertEquals( 4, s.getNumeratorDof().intValue() );
833         assertEquals( 25.14, s.getF(), 0.01 );
834         assertEquals( 1.743e-05, s.getP(), 0.001 );
835         anova = s.getAnova();
836         assertEquals( 38.411, anova.getMainEffectF( "Treatment" ), 0.001 );
837 
838         s = sums.get( "1456759_at" );
839         assertNotNull( s.getContrastCoefficients() );
840         assertEquals( 11, s.getResidualDof().intValue() );
841         assertEquals( 4, s.getNumeratorDof().intValue() );
842         assertEquals( 7.903, s.getF(), 0.01 );
843         assertEquals( 0.002960, s.getP(), 0.001 );
844         anova = s.getAnova();
845         assertEquals( 10.3792, anova.getMainEffectF( "Treatment" ), 0.001 );
846         assertEquals( 2.6253, anova.getMainEffectF( "Genotype" ), 0.001 );
847     }
848 
849     /**
850      * Sanity check.
851      *
852      * @throws Exception
853      */
854     @Test
855     public void testTwoWayAnovaUnfittable() throws Exception {
856         DoubleMatrixReader f = new DoubleMatrixReader();
857         DoubleMatrix<String, String> testMatrix = f.read( this.getClass()
858                 .getResourceAsStream( "/data/lmtest10.dat.txt" ) );
859 
860         StringMatrixReader of = new StringMatrixReader();
861         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
862                 "/data/lmtest10.des.txt" ) );
863 
864         DesignMatrix d = new DesignMatrix( sampleInfo, true );
865 
866         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix );
867         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
868         assertEquals( 1, sums.size() );
869 
870         for ( LinearModelSummary lms : sums.values() ) {
871             GenericAnovaResult a = lms.getAnova();
872             assertNotNull( a );
873             assertEquals( Double.NaN, a.getMainEffectP( "CellType" ), 0.0001 );
874             assertEquals( Double.NaN, a.getMainEffectP( "SamplingTimePoint" ), 0.0001 );
875         }
876 
877         DoubleMatrix2D coefficients = fit.getCoefficients();
878         DoubleMatrix2D residuals = fit.getResiduals();
879 
880         assertEquals( 2.238, coefficients.get( 0, 0 ), 0.0001 ); // mean.
881         assertEquals( 0.0, coefficients.get( 1, 0 ), 0.0001 );
882 
883         for ( int i = 0; i < residuals.rows(); i++ ) {
884             assertEquals( 0.0, residuals.get( 0, i ), 0.00001 );
885         }
886     }
887 
888     /**
889      * Check for problem reported by TF -- Gemma gives slightly different result. Problem is not at this level.
890      *
891      * @throws Exception
892      */
893     @Test
894     public void testTwoWayAnovaWithInteractions() throws Exception {
895         DoubleMatrixReader f = new DoubleMatrixReader();
896         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
897                 "/data/GSE8441_expmat_8probes.txt" ) );
898 
899         StringMatrixReader of = new StringMatrixReader();
900         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
901                 "/data/606_GSE8441_expdesign.data.txt" ) );
902         DesignMatrix d = new DesignMatrix( sampleInfo, true );
903         d.addInteraction();
904 
905         assertEquals( 4, d.getMatrix().columns() );
906 
907         assertEquals( 22, testMatrix.columns() );
908 
909         LeastSquaresFit fit = new LeastSquaresFit( d, testMatrix );
910 
911         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
912         assertEquals( 8, sums.size() );
913 
914         for ( LinearModelSummary lms : sums.values() ) {
915             GenericAnovaResult a = lms.getAnova();
916             assertNotNull( a );
917         }
918 
919         LinearModelSummary s = sums.get( "217757_at" );
920         GenericAnovaResult anova = s.getAnova();
921 
922         assertNotNull( s.getContrastCoefficients() );
923         assertEquals( 0.763, s.getContrastCoefficients().get( 2, 3 ), 0.001 );
924         assertEquals( 18, s.getResidualDof().intValue() );
925         assertEquals( 3, s.getNumeratorDof().intValue() );
926         assertEquals( 0.299, s.getF(), 0.01 );
927         assertEquals( 0.8257, s.getP(), 0.001 );
928 
929         assertEquals( 0.5876, anova.getMainEffectF( "Treatment" ), 0.0001 );
930         assertEquals( 0.5925, anova.getInteractionEffectP(), 0.001 );
931 
932         s = sums.get( "202851_at" );
933         anova = s.getAnova();
934         assertNotNull( s.getContrastCoefficients() );
935         assertEquals( 0.787, s.getContrastCoefficients().get( 2, 3 ), 0.001 );
936         assertEquals( 18, s.getResidualDof().intValue() );
937         assertEquals( 3, s.getNumeratorDof().intValue() );
938         assertEquals( 0.1773, s.getF(), 0.01 );
939         assertEquals( 0.9104, s.getP(), 0.001 );
940 
941         assertEquals( 0.3777, anova.getMainEffectF( "Treatment" ), 0.0001 );
942         assertEquals( 0.9956, anova.getInteractionEffectP(), 0.001 );
943     }
944 
945     /**
946      * No missing values; Two-way ANOVA with interaction, PLUS a continuous covariate.
947      *
948      * @throws Exception
949      */
950     @Test
951     public void testTwoWayTwoLevelsOneContinousInteractionC() throws Exception {
952         DoubleMatrixReader f = new DoubleMatrixReader();
953         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
954                 "/data/example.madata.small.txt" ) );
955 
956         ObjectMatrix<String, String, Object> design = new ObjectMatrixImpl<>( 9, 3 );
957 
958         design.set( 0, 0, "A" );
959         design.set( 1, 0, "A" );
960         design.set( 2, 0, "A" );
961         design.set( 3, 0, "A" );
962         design.set( 4, 0, "B" );
963         design.set( 5, 0, "B" );
964         design.set( 6, 0, "B" );
965         design.set( 7, 0, "B" );
966         design.set( 8, 0, "B" );
967         design.set( 0, 1, 0.12 );
968         design.set( 1, 1, 0.24 );
969         design.set( 2, 1, 0.48 );
970         design.set( 3, 1, 0.96 );
971         design.set( 4, 1, 0.12 );
972         design.set( 5, 1, 0.24 );
973         design.set( 6, 1, 0.48 );
974         design.set( 7, 1, 0.96 );
975         design.set( 8, 1, 0.96 );
976         design.set( 0, 2, "C" );
977         design.set( 1, 2, "C" );
978         design.set( 2, 2, "D" );
979         design.set( 3, 2, "D" );
980         design.set( 4, 2, "C" );
981         design.set( 5, 2, "C" );
982         design.set( 6, 2, "D" );
983         design.set( 7, 2, "D" );
984         design.set( 8, 2, "D" );
985         design.addColumnName( "Treat" );
986         design.addColumnName( "Value" );
987         design.addColumnName( "Geno" );
988 
989         DesignMatrix designMatrix = new DesignMatrix( design, true );
990         designMatrix.addInteraction( "Treat", "Geno" );
991         LeastSquaresFit fit = new LeastSquaresFit( designMatrix, testMatrix );
992 
993         Map<String, LinearModelSummary> sums = fit.summarizeByKeys( true );
994 
995         LinearModelSummary s = sums.get( "1553129_at" );
996         assertEquals( 1.791, s.getF(), 0.01 );
997         assertEquals( 0.2930, s.getP(), 0.001 );
998         assertEquals( 3.71542, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
999         assertEquals( 0.184, s.getContrastCoefficients().get( 1, 3 ), 0.001 ); // this ordering might change?
1000         assertEquals( 0.137, s.getContrastCoefficients().get( 4, 3 ), 0.001 ); // this ordering might change?
1001         GenericAnovaResult a = s.getAnova();
1002         assertEquals( 0.137, a.getInteractionEffectP(), 0.001 );
1003 
1004         s = sums.get( "232018_at" );
1005         assertEquals( 0.7167, s.getF(), 0.01 );
1006         assertEquals( 0.6235, s.getP(), 0.001 );
1007         assertEquals( 6.8873, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
1008         assertEquals( 0.587932, s.getContrastCoefficients().get( 1, 3 ), 0.001 );
1009         a = s.getAnova();
1010         assertEquals( 0.2904, a.getInteractionEffectP(), 0.001 );
1011 
1012         for ( LinearModelSummary lms : sums.values() ) {
1013             GenericAnovaResult anova = lms.getAnova();
1014             assertNotNull( anova );
1015             Double interactionEffectP = anova.getInteractionEffectP();
1016             assertNotNull( interactionEffectP );
1017             assertTrue( !Double.isNaN( interactionEffectP ) );
1018 
1019         }
1020 
1021     }
1022 
1023     @Test
1024     public void testVectorRegress() {
1025 
1026         DoubleMatrix1D vectorA = new DenseDoubleMatrix1D( new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 } );
1027         DoubleMatrix1D vectorB = new DenseDoubleMatrix1D( new double[] { 1, 2, 2, 3, 3, 4, 4, 5, 5, 6 } );
1028 
1029         LeastSquaresFit fit = new LeastSquaresFit( vectorA, vectorB );
1030 
1031         DoubleMatrix2D coefficients = fit.getCoefficients();
1032         DoubleMatrix2D residuals = fit.getResiduals();
1033 
1034         assertEquals( 0.666666, coefficients.get( 0, 0 ), 0.0001 );
1035         assertEquals( 0.5152, coefficients.get( 1, 0 ), 0.0001 );
1036 
1037         double[] expectedResiduals = new double[] { -0.1818182, 0.3030303, -0.2121212, 0.2727273, -0.2424242,
1038                 0.2424242, -0.2727273, 0.2121212, -0.3030303, 0.1818182 };
1039 
1040         for ( int i = 0; i < expectedResiduals.length; i++ ) {
1041             assertEquals( expectedResiduals[i], residuals.get( 0, i ), 0.00001 );
1042         }
1043 
1044     }
1045 
1046     /**
1047      * @throws Exception
1048      */
1049     @Test
1050     public void testVectorWeightedRegress() throws Exception {
1051         // a<-c( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 )
1052         //  b<-c(1, 2, 2, 3, 3, 4, 4, 5, 5, 6)
1053         //  w<-1/a
1054         DoubleMatrix1D vectorA = new DenseDoubleMatrix1D( new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 } );
1055         DoubleMatrix1D vectorB = new DenseDoubleMatrix1D( new double[] { 1, 2, 2, 3, 3, 4, 4, 5, 5, 6 } );
1056         DoubleMatrix1D w = vectorA.copy().assign( Functions.inv );
1057 
1058         LeastSquaresFit fit = new LeastSquaresFit( vectorA, vectorB, w );
1059 
1060         DoubleMatrix2D coefficients = fit.getCoefficients();
1061         DoubleMatrix2D residuals = fit.getResiduals();
1062         //lm(b~a, weights=w)
1063         //lm.wfit(model.matrix(~a, factor(a)), b, w)
1064         assertEquals( 0.60469, coefficients.get( 0, 0 ), 0.0001 );
1065         assertEquals( 0.52642, coefficients.get( 1, 0 ), 0.0001 );
1066 
1067         double[] expectedResiduals = new double[] { -0.1311097, 0.3424702, -0.1839499, 0.2896301, -0.2367900,
1068                 0.2367900, -0.2896301, 0.1839499, -0.3424702, 0.1311097 };
1069 
1070         for ( int i = 0; i < expectedResiduals.length; i++ ) {
1071             assertEquals( expectedResiduals[i], residuals.get( 0, i ), 0.00001 );
1072         }
1073 
1074         double[] expectedFitted = new double[] { 1.13111, 1.65753, 2.18395, 2.71037, 3.23679, 3.76321, 4.28963,
1075                 4.81605, 5.34247, 5.86889 };
1076 
1077         for ( int i = 0; i < expectedFitted.length; i++ ) {
1078             assertEquals( expectedFitted[i], fit.getFitted().get( 0, i ), 0.00001 );
1079         }
1080     }
1081 
1082     /**
1083      * @throws Exception
1084      */
1085     @Test
1086     public void testVectorWeightedRegressWithMissing() throws Exception {
1087 
1088         DoubleMatrixReader f = new DoubleMatrixReader();
1089         DoubleMatrix<String, String> testMatrix = f.read( this.getClass().getResourceAsStream(
1090                 "/data/example.madata.withmissing.small.txt" ) );
1091         DoubleMatrix1D libSize = MatrixStats.colSums( testMatrix );
1092         testMatrix = MatrixStats.convertToLog2Cpm( testMatrix, libSize );
1093 
1094         StringMatrixReader of = new StringMatrixReader();
1095         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
1096                 "/data/example.metadata.small.txt" ) );
1097         DesignMatrix designMatrix = new DesignMatrix( sampleInfo );
1098         DoubleMatrix2D weights = new DenseDoubleMatrix2D( testMatrix.asArray() );
1099         weights.assign( Functions.inv );
1100 
1101         LeastSquaresFit fit = new LeastSquaresFit( designMatrix, testMatrix, weights );
1102         assertTrue( fit.isHasMissing() );
1103 
1104         DoubleMatrix2D coefficients = fit.getCoefficients();
1105         DoubleMatrix2D residuals = fit.getResiduals();
1106 
1107         assertEquals( 15.339801, coefficients.get( 0, 0 ), 0.0001 );
1108         assertEquals( -0.024058, coefficients.get( 1, 1 ), 0.0001 );
1109         assertEquals( -0.059586, coefficients.get( 2, 18 ), 0.0001 );
1110 
1111         assertEquals( -0.073732, residuals.get( 0, 0 ), 0.0001 );
1112         assertEquals( -0.064656, residuals.get( 1, 1 ), 0.0001 );
1113         assertEquals( -0.085214, residuals.get( 18, 8 ), 0.0001 );
1114         assertTrue( Double.isNaN( residuals.get( 4, 2 ) ) );
1115 
1116     }
1117 
1118     /**
1119      * Tests limma-like functionality
1120      * 
1121      * Multiple levels per factor, unbalanced design
1122      * 
1123      * @throws Exception
1124      */
1125     @Test
1126     public void testNHBE() throws Exception {
1127         DoubleMatrixReader f = new DoubleMatrixReader();
1128         DoubleMatrix<String, String> testMatrix = f.read( new GZIPInputStream( this.getClass().getResourceAsStream(
1129                 "/data/NHBE_transcriptome_data.txt.gz" ) ) );
1130 
1131         StringMatrixReader of = new StringMatrixReader();
1132         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
1133                 "/data/NHBE_design.txt" ) );
1134 
1135         DesignMatrix designMatrix = new DesignMatrix( sampleInfo );
1136         designMatrix.addInteraction();
1137         designMatrix.setBaseline( "time", "1_h" );
1138         designMatrix.setBaseline( "Treatment", "control" );
1139 
1140         LeastSquaresFit fit = new LeastSquaresFit( designMatrix, testMatrix );
1141 
1142         //  System.err.println( designMatrix );
1143         //   List<LinearModelSummary> sums = fit.summarize( true );
1144 
1145         ModeratedTstat.ebayes( fit );
1146 
1147         /////////////
1148         //  System.err.println( "------- After ebayes ------" );
1149         List<LinearModelSummary> sums = fit.summarize( true );
1150 
1151         // fit3$sigma[1]
1152         assertEquals( 0.34927, sums.get( 0 ).getSigma(), 0.0001 );
1153         assertEquals( 1.3859, sums.get( 0 ).getPriorDof(), 0.01 );
1154 
1155         /*
1156          * TODO: add more tests here.
1157          */
1158     }
1159 
1160     /**
1161      * Tests limma-like functionality. Balanced 2x2 design
1162      * 
1163      * @throws Exception
1164      */
1165     @Test
1166     public void testEstrogen() throws Exception {
1167         DoubleMatrixReader f = new DoubleMatrixReader();
1168         DoubleMatrix<String, String> testMatrix = f.read( new GZIPInputStream( this.getClass().getResourceAsStream(
1169                 "/data/estrogen.data.txt.gz" ) ) );
1170 
1171         StringMatrixReader of = new StringMatrixReader();
1172         StringMatrix<String, String> sampleInfo = of.read( this.getClass().getResourceAsStream(
1173                 "/data/estrogen.meta.txt" ) );
1174 
1175         DesignMatrix designMatrix = new DesignMatrix( sampleInfo );
1176         designMatrix.addInteraction();
1177         LeastSquaresFit fit = new LeastSquaresFit( designMatrix, testMatrix );
1178 
1179         //  System.err.println( designMatrix );
1180         List<LinearModelSummary> sums = fit.summarize( true );
1181         //  System.err.println( fit.getCoefficients().viewColumn( 0 ) );
1182 
1183         LinearModelSummary s = sums.get( 0 );
1184         assertEquals( 3.8976, s.getF(), 0.01 );
1185         assertEquals( 0.11092, s.getP(), 0.001 );
1186         assertEquals( 9.69220, s.getContrastCoefficients().get( 0, 0 ), 0.001 );
1187         assertEquals( -1.4517, s.getContrastCoefficients().get( 1, 2 ), 0.001 ); // tstat
1188         assertEquals( 0.220, s.getContrastCoefficients().get( 1, 3 ), 0.001 ); // pvalue
1189         assertEquals( 4, s.getResidualDof().intValue() );
1190         GenericAnovaResult a = s.getAnova();
1191         assertEquals( 0.24143, a.getMainEffectF( "time" ), 0.0001 );
1192         assertEquals( 1, a.getMainEffectDof( "time" ).intValue() );
1193         assertEquals( 4, a.getResidualDf().intValue() );
1194         assertEquals( 2.43873, a.getInteractionEffectF(), 0.001 );
1195         assertEquals( 0.19340, a.getInteractionEffectP(), 0.001 );
1196 
1197         /////////////
1198         ModeratedTstat.ebayes( fit );
1199         sums = fit.summarize( true );
1200         LinearModelSummary x = sums.get( 0 );
1201 
1202         assertEquals( 0.0765, x.getSigma(), 0.0001 );
1203         assertEquals( 4.48, x.getPriorDof(), 0.01 ); // we get 4.479999 or something.
1204 
1205         assertEquals( 3.8976, x.getF(), 0.01 );
1206         assertEquals( 0.11092, x.getP(), 0.001 );
1207         // topTable(fit3, coef=2,n=dim(fit3)[1], sort.by="none" )["100_g_at",]
1208         assertEquals( 9.692196, x.getContrastCoefficients().get( 0, 0 ), 0.0001 );
1209         assertEquals( -0.92608, x.getContrastCoefficients().get( 1, 2 ), 0.001 ); // tstat
1210         assertEquals( 0.38, x.getContrastCoefficients().get( 1, 3 ), 0.001 ); // pvalue
1211         assertEquals( 4, x.getResidualDof().intValue() );
1212         // topTable(fit3, coef=4,n=dim(fit3)[1], sort.by="none" )["100_g_at",]
1213         assertEquals( 0.34671, x.getContrastCoefficients().get( 3, 3 ), 0.0001 ); // interaction pvalue
1214 
1215         GenericAnovaResult ax = x.getAnova();
1216         // classifyTestsF(fit3$t[,c(2)], df=fit3$df.residual, cor.matrix=cov2cor(fit3$cov.coefficients[2,2]), fstat.only = T)[1]
1217         assertEquals( 0.098252, ax.getMainEffectF( "time" ), 0.0001 );
1218         assertEquals( 1, ax.getMainEffectDof( "time" ).intValue() );
1219         assertEquals( 8.48, ax.getResidualDf(), 0.001 );
1220 
1221         assertEquals( 3.6678, ax.getMainEffectF( "dose" ), 0.0001 );
1222         assertEquals( 1, ax.getMainEffectDof( "dose" ).intValue() );
1223         assertEquals( 8.48, ax.getResidualDf(), 0.001 ); // 4 + 4.48
1224 
1225         // sum(f4^2)/1/sqig^2
1226         // pf(0.99247, 1, 8.48, lower.tail=F) 
1227         assertEquals( 0.99247, ax.getInteractionEffectF(), 0.0001 );
1228         assertEquals( 0.34671, ax.getInteractionEffectP(), 0.0001 );
1229 
1230     }
1231 }