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