Source code

001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.survival.cox;
022
023import org.biojava.nbio.survival.cox.stats.ChiSq;
024import org.biojava.nbio.survival.kaplanmeier.figure.ExpressionFigure;
025import org.biojava.nbio.survival.kaplanmeier.figure.KaplanMeierFigure;
026
027import java.text.DecimalFormat;
028import java.util.ArrayList;
029import java.util.LinkedHashMap;
030
031//import org.biojava.nbio.survival.cox.comparators.SurvivalInfoComparator;
032
033/**
034 * Holds the results of a cox analysis where calling dump(), toString() will give an output similar to R
035 * @author Scooter Willis <willishf at gmail dot com>
036 */
037public class CoxInfo {
038
039        private WaldTestInfo waldTestInfo = null;
040        String message = "";
041        Integer maxIterations = null;
042        Double eps = null;
043        Double toler = null;
044        CoxMethod method;
045        private double[][] imat = null; //:the variance matrix at beta=final
046        private double[][] naive_imat = null; //the original variance matrix used in residuals
047        double[] u = new double[0]; //:score vector
048        int iterations = 0; //:actual number of iterations used
049        int flag = 0; // success flag  1000  did not converge 1 to nvar: rank of the solution
050        double logTest = 0;
051        double logTestpval = 0;
052        double loglikInit = 0; //loglik at beta=initial values, at beta=final
053        double loglikFinal = 0;
054        Double scoreLogrankTest;
055        private Double rscore = null; //robust score
056        private Double rscoreLogrankTestpvalue = null;
057        private double degreeFreedom;
058        private Double scoreLogrankTestpvalue;
059        int numSamples = 0;
060        int numEvents = 0;
061        private LinkedHashMap<String, String> metaDataFilter = null;
062        private LinkedHashMap<String, CoxCoefficient> coefficientsList = new LinkedHashMap<String, CoxCoefficient>();
063        LinkedHashMap<Double, Double> baselineSurvivorFunction = new LinkedHashMap<Double, Double>();
064        ArrayList<SurvivalInfo> survivalInfoList = new ArrayList<SurvivalInfo>();
065        /**
066         *
067         */
068        public KaplanMeierFigure kmf = null;
069        /**
070         *
071         */
072        public ExpressionFigure ef = null;
073
074        /**
075         *
076         * @return
077         */
078        public ArrayList<SurvivalInfo> getSurvivalInfoList() {
079                return survivalInfoList;
080        }
081
082        /**
083         *
084         * @param var
085         * @throws Exception
086         */
087        public void setVariance(double[][] var) throws Exception {
088
089                //   if (Math.abs(var[0][1] - var[1][0]) > .0000000000001) { //in the CoxCC correction looks like a precision error keeps these from being equal 10-19
090                //       throw new Exception("Expecting diagonal to be equal");
091                //   }
092                imat = new double[var.length][var[0].length];
093                for (int i = 0; i < var.length; i++) {
094                        for (int j = 0; j < var[0].length; j++) {
095                                imat[i][j] = var[i][j];
096                        }
097                }
098                calcSummaryValues();
099        }
100
101        /**
102         *
103         * @return
104         */
105        public double[][] getVariance() {
106                double[][] var = new double[imat.length][imat[0].length];
107                for (int i = 0; i < var.length; i++) {
108                        for (int j = 0; j < var[0].length; j++) {
109                                var[i][j] = imat[i][j];
110                        }
111                }
112
113                return var;
114        }
115
116        /**
117         *
118         * @param var
119         * @throws Exception
120         */
121        public void setNaiveVariance(double[][] var) throws Exception {
122//        if (var[0][1] != var[1][0]) {
123//            throw new Exception("Expecting diagonal to be equal");
124//        }
125
126
127                naive_imat = new double[var.length][var[0].length];
128                for (int i = 0; i < var.length; i++) {
129                        for (int j = 0; j < var[0].length; j++) {
130                                naive_imat[i][j] = var[i][j];
131                        }
132                }
133
134                calcSummaryValues();
135        }
136
137        /**
138         *
139         * @return
140         */
141        public double[][] getNaiveVariance() {
142                double[][] var = new double[imat.length][imat[0].length];
143                for (int i = 0; i < var.length; i++) {
144                        for (int j = 0; j < var[0].length; j++) {
145                                var[i][j] = naive_imat[i][j];
146                        }
147                }
148
149                return var;
150
151        }
152
153        /**
154         *
155         * @param data
156         */
157        public void setSurvivalInfoList(ArrayList<SurvivalInfo> data) {
158                survivalInfoList = data;
159                numSamples = data.size();
160
161                for (SurvivalInfo si : data) {
162                        if (si.getStatus() == 1) {
163                                numEvents++;
164                        }
165                }
166        }
167
168        /**
169         *
170         * @return
171         */
172        public double[] getWeighted() {
173                double[] weighted = new double[survivalInfoList.size()];
174                int p = 0;
175                for (SurvivalInfo si : this.survivalInfoList) {
176                        weighted[p] = si.getWeight();
177                        p++;
178                }
179                return weighted;
180        }
181
182        /**
183         *
184         * @return
185         */
186        public double[][] getVariableResiduals() {
187                ArrayList<String> variables = new ArrayList<String>(coefficientsList.keySet());
188                double[][] rr = new double[survivalInfoList.size()][variables.size()];
189                int p = 0;
190                for (SurvivalInfo si : this.survivalInfoList) {
191                        int i = 0;
192                        for (String v : variables) {
193                                rr[p][i] = si.getResidualVariable(v);
194                                i++;
195                        }
196                        p++;
197                }
198
199                return rr;
200        }
201
202        /**
203         *
204         * @param rr
205         */
206        public void setVariableResiduals(double[][] rr) {
207                ArrayList<String> variables = new ArrayList<String>(coefficientsList.keySet());
208
209                int p = 0;
210                for (SurvivalInfo si : this.survivalInfoList) {
211                        int i = 0;
212                        for (String v : variables) {
213                                si.setResidualVariable(v, rr[p][i]);
214                                i++;
215                        }
216                        p++;
217                }
218
219        }
220
221        /**
222         *
223         * @return
224         */
225        public int getNumberCoefficients() {
226                return coefficientsList.size();
227        }
228
229        /**
230         *
231         * @param name
232         * @return
233         */
234        public CoxCoefficient getCoefficient(String name) {
235                return coefficientsList.get(name);
236        }
237
238        /**
239         *
240         * @param name
241         * @param coefficient
242         */
243        public void setCoefficient(String name, CoxCoefficient coefficient) {
244                coefficientsList.put(name, coefficient);
245        }
246
247        /**
248         *
249         * @param header
250         * @param beginLine
251         * @param beginCell
252         * @param endCell
253         * @param endLine
254         * @return
255         */
256        public String getCoefficientText(boolean header, String beginLine, String beginCell, String endCell, String endLine) {
257                String o = "";
258                if (header) {
259                        String robust = "";
260                        if (naive_imat != null) {
261                                robust = beginCell + "robust se" + endCell;
262                        }
263                        o = o + beginLine + beginCell + fmtpl("", 9) + endCell + beginCell + fmtpl("coef", 9) + endCell + beginCell + fmtpl("se(coef)", 9) + endCell + robust + beginCell + fmtpl("z", 9) + endCell + beginCell + fmtpl("p-value", 9) + endCell + beginCell + fmtpl("HR", 9) + endCell + beginCell + fmtpl("lower .95", 9) + endCell + beginCell + fmtpl("upper .95", 9) + endCell + endLine;
264                }//Coefficients,Coe,StdErr,HR,p-value,HR Lo 95%,HR Hi 95%
265
266                for (CoxCoefficient coe : coefficientsList.values()) {
267                        String robust = "";
268                        String stderror = "";
269                        if (naive_imat != null) {
270                                stderror = beginCell + fmt(coe.getRobustStdError(), 5, 9) + endCell;
271                                robust = beginCell + fmt(coe.getStdError(), 5, 9) + endCell;
272                        } else {
273                                stderror = beginCell + fmt(coe.getStdError(), 5, 9) + endCell;
274                        }
275                        o = o + beginLine + beginCell + fmtpr(coe.getName(), 9) + endCell + beginCell + fmt(coe.getCoeff(), 5, 9) + stderror + robust + endCell + beginCell + fmt(coe.getZ(), 5, 9) + endCell + beginCell + fmt(coe.getPvalue(), 6, 9) + endCell + beginCell + fmt(coe.getHazardRatio(), 3, 9) + endCell + beginCell + fmt(coe.getHazardRatioLoCI(), 3, 9) + endCell + beginCell + fmt(coe.getHazardRatioHiCI(), 3, 9) + endCell + endLine;
276                }
277                return o;
278        }
279
280        /**
281         *
282         * @param d
283         * @param precision
284         * @param pad
285         * @return
286         */
287        public static String fmt(Double d, int precision, int pad) {
288                if(d == null)
289                        return "";
290                if(Double.isNaN(d))
291                        return "";
292                String value = "";
293                DecimalFormat dfe = new DecimalFormat("0.00E0");
294                String dpad = "0.";
295                double p = 1.0;
296                for (int i = 0; i < (precision); i++) {
297                        dpad = dpad + "0";
298                        p = p / 10.0;
299                }
300                DecimalFormat df = new DecimalFormat(dpad);
301                if (Math.abs(d) >= p) {
302                        value = df.format(d);
303                } else {
304                        value = dfe.format(d);
305                }
306                int length = value.length();
307                int extra = pad - length;
308                if (extra > 0) {
309                        for (int i = 0; i < extra; i++) {
310                                value = " " + value;
311                        }
312                }
313                return value;
314        }
315
316        /**
317         *
318         */
319        private void calcSummaryValues() {
320
321                //beta
322
323                ArrayList<String> variables = new ArrayList<String>(coefficientsList.keySet());
324                for (int i = 0; i < variables.size(); i++) {
325                        String variable = variables.get(i);
326                        CoxCoefficient coe = coefficientsList.get(variable);
327                        coe.setStdError(Math.sqrt(imat[i][i])); //values can be updated to reflect new error
328                        if (naive_imat != null) {
329                                coe.setRobustStdError(Math.sqrt(naive_imat[i][i]));
330                        }
331                        coe.setZ(coe.getCoeff() / coe.getStdError());
332                        coe.setPvalue(ChiSq.norm(Math.abs(coe.getCoeff() / coe.getStdError())));
333                        //z <- qnorm((1 + conf.int)/2, 0, 1)
334                        double z = 1.959964;
335                        coe.setHazardRatioLoCI(Math.exp(coe.getCoeff() - z * coe.getStdError()));
336                        coe.setHazardRatioHiCI(Math.exp(coe.getCoeff() + z * coe.getStdError()));
337                }
338
339                logTest = -2 * (loglikInit - loglikFinal);
340                logTestpval = ChiSq.chiSq(logTest, (int) degreeFreedom);
341
342                scoreLogrankTestpvalue = ChiSq.chiSq(scoreLogrankTest, (int) degreeFreedom);
343                if (rscore != null) {
344                        rscoreLogrankTestpvalue = ChiSq.chiSq(rscore, (int) degreeFreedom);
345                }
346        }
347
348        /**
349         *
350         */
351        public void dump() {
352
353                //need an ordered list for comparing to R dumps
354
355//        ArrayList<SurvivalInfo> orderedSurvivalInfoList = new ArrayList<SurvivalInfo>(survivalInfoList);
356//        SurvivalInfoComparator sicSort = new SurvivalInfoComparator();
357                //       Collections.sort(orderedSurvivalInfoList,sicSort);
358
359
360                System.out.println();
361                System.out.println("$coef");
362                for (CoxCoefficient coe : coefficientsList.values()) {
363                        System.out.print(coe.getCoeff() + " ");
364                }
365                System.out.println();
366                System.out.println("$means");
367
368                for (CoxCoefficient coe : coefficientsList.values()) {
369                        System.out.print(coe.getMean() + " ");
370                }
371                System.out.println();
372                System.out.println("$u");
373
374                for (double d : u) {
375                        System.out.print(d + " ");
376                }
377
378                System.out.println();
379                System.out.println("$imat");
380                for (int i = 0; i < imat.length; i++) {
381                        for (int j = 0; j < imat[0].length; j++) {
382                                System.out.print(imat[i][j] + " ");
383                        }
384                        System.out.println();
385                }
386
387                if (this.naive_imat != null) {
388                        System.out.println("$naive_imat");
389                        for (int i = 0; i < naive_imat.length; i++) {
390                                for (int j = 0; j < naive_imat[0].length; j++) {
391                                        System.out.print(naive_imat[i][j] + " ");
392                                }
393                                System.out.println();
394                        }
395                }
396
397                System.out.println();
398                System.out.println("$loglik");
399
400                System.out.println(loglikInit + " " + loglikFinal);
401
402                System.out.println();
403                System.out.println("$sctest");
404
405                System.out.println(this.scoreLogrankTest);
406
407                System.out.println("$iter");
408                System.out.println(this.iterations);
409
410                System.out.println("$flag");
411                System.out.println(flag);
412
413                System.out.println();
414//        if (false) {
415//            System.out.println("ID      LP       Score      Residuals");
416//            for (SurvivalInfo si : orderedSurvivalInfoList) {
417//                System.out.println(si.getOrder() + " " + si.getLinearPredictor() + " " + si.getScore() + " " + si.getResidual());
418//
419//            }
420//            System.out.println();
421//            ArrayList<String> variables = new ArrayList<String>(coefficientsList.keySet());
422//            System.out.print("Sample");
423//            for (String v : variables) {
424//                System.out.print("    " + v);
425//            }
426//            System.out.println("rr");
427//            for (SurvivalInfo si : orderedSurvivalInfoList) {
428//                System.out.print(si.getOrder());
429//                for (String v : variables) {
430//                    System.out.print("   " + si.getResidualVariable(v));
431//                }
432//                System.out.println();
433//            }
434//        }
435
436        }
437
438        /**
439         *
440         * @param d
441         * @param pad
442         * @return
443         */
444        public String fmtpr(String d, int pad) {
445                int length = d.length();
446                int extra = pad - length;
447                if (extra < 0) {
448                        extra = 0;
449                }
450                String v = d;
451                for (int i = 0; i < extra; i++) {
452                        v = v + " ";
453                }
454                return v;
455        }
456
457        /**
458         * Pad left a string with spaces
459         *
460         * @param d
461         * @param pad
462         * @return
463         */
464        public String fmtpl(String d, int pad) {
465                int length = d.length();
466                int extra = pad - length;
467                if (extra < 0) {
468                        extra = 0;
469                }
470                String v = d;
471                for (int i = 0; i < extra; i++) {
472                        v = " " + v;
473                }
474                return v;
475        }
476
477        @Override
478        public String toString() {
479                return toString("", " ", "\r\n");
480        }
481
482        /**
483         *
484         * @param beginLine
485         * @param del
486         * @param endLine
487         * @return
488         */
489        public String toString(String beginLine, String del, String endLine) {
490
491
492
493                String o = beginLine + fmtpl("", 9) + fmtpl("Avg", 9) + fmtpl("SD", 9) + endLine;
494                for (CoxCoefficient coe : coefficientsList.values()) {
495                        o = o + beginLine + fmtpr(coe.getName(), 9) + fmt(coe.getMean(), 4, 9) + fmt(coe.getStandardDeviation(), 4, 9) + endLine;
496
497                }
498
499                o = o + beginLine + endLine;
500
501
502
503                o = o + beginLine + "n= " + this.numSamples + ", number of events=" + this.numEvents + endLine;
504                o = o + getCoefficientText(true, beginLine, del, "", endLine);
505
506                o = o + beginLine + endLine;
507
508                if (baselineSurvivorFunction.size() > 0) {
509                        o = o + beginLine + "Baseline Survivor Function (at predictor means)" + endLine;
510                        for (Double time : baselineSurvivorFunction.keySet()) {
511                                Double mean = baselineSurvivorFunction.get(time);
512                                o = o + beginLine + fmt(time, 4, 10) + fmt(mean, 4, 10) + endLine;
513                        }
514                }
515
516                o = o + beginLine + endLine;
517                o = o + beginLine + "Overall Model Fit" + endLine;
518                o = o + beginLine + "Iterations=" + iterations + endLine;
519
520
521                o = o + beginLine + "Likelihood ratio test = " + fmt(this.logTest, 2, 0) + " df=" + this.degreeFreedom + " p-value=" + fmt(this.logTestpval, 7, 0) + endLine;
522
523                o = o + beginLine + "Wald test             = " + fmt(waldTestInfo.getTest(), 2, 0) + " df=" + waldTestInfo.getDf() + " p-value=" + fmt(waldTestInfo.getPvalue(), 7, 0) + endLine;
524                o = o + beginLine + "Score (logrank) test  = " + fmt(scoreLogrankTest, 2, 0) + " df=" + ((int) (this.degreeFreedom)) + " p-value=" + fmt(this.scoreLogrankTestpvalue, 7, 0);
525
526                if (this.rscore != null) {
527                        o = o + ",  Robust = " + fmt(rscore, 2, 0) + " p-value=" + fmt(rscoreLogrankTestpvalue, 7, 0);
528
529                }
530
531                o = o + endLine;
532
533                //       o = o + "Rank of solution flag=" + flag + "\r\n";
534                //       o = o + "Log lik Initial=" + loglikInit + "\r\n";
535                //       o = o + "Log lik Final=" + loglikFinal + "\r\n";
536                o = o + beginLine + "Method=" + method.name() + endLine;
537
538
539
540
541                return o;
542        }
543
544        /**
545         * @return the scoreLogrankTest
546         */
547        public double getChiSquare() {
548                return scoreLogrankTest;
549        }
550
551        /**
552         * @return the degreeFreedom
553         */
554        public double getDegreeFreedom() {
555                return degreeFreedom;
556        }
557
558        /**
559         * @return the scoreLogrankTestpvalue
560         */
561        public double getOverallModelFitPvalue() {
562                return scoreLogrankTestpvalue;
563        }
564
565        /**
566         * @return the rscore
567         */
568        public Double getRscore() {
569                return rscore;
570        }
571
572        /**
573         * @param rscore the rscore to set
574         */
575        public void setRscore(Double rscore) {
576                this.rscore = rscore;
577                if (rscore != null) {
578                        rscoreLogrankTestpvalue = ChiSq.chiSq(rscore, (int) degreeFreedom);
579                }
580        }
581
582        /**
583         * @return the rscoreLogrankTestpvalue
584         */
585        public Double getRscoreLogrankTestpvalue() {
586                return rscoreLogrankTestpvalue;
587        }
588
589        /**
590         * @param rscoreLogrankTestpvalue the rscoreLogrankTestpvalue to set
591         */
592        public void setRscoreLogrankTestpvalue(Double rscoreLogrankTestpvalue) {
593                this.rscoreLogrankTestpvalue = rscoreLogrankTestpvalue;
594        }
595
596        /**
597         * @return the scoreLogrankTest
598         */
599        public Double getScoreLogrankTest() {
600                return scoreLogrankTest;
601        }
602
603        /**
604         * @param scoreLogrankTest the scoreLogrankTest to set
605         */
606        public void setScoreLogrankTest(Double scoreLogrankTest) {
607                this.scoreLogrankTest = scoreLogrankTest;
608        }
609
610        /**
611         * @return the scoreLogrankTestpvalue
612         */
613        public Double getScoreLogrankTestpvalue() {
614                return scoreLogrankTestpvalue;
615        }
616
617        /**
618         * @param scoreLogrankTestpvalue the scoreLogrankTestpvalue to set
619         */
620        public void setScoreLogrankTestpvalue(Double scoreLogrankTestpvalue) {
621                this.scoreLogrankTestpvalue = scoreLogrankTestpvalue;
622        }
623
624        /**
625         * @return the metaDataFilter
626         */
627        public LinkedHashMap<String, String> getMetaDataFilter() {
628                return metaDataFilter;
629        }
630
631        /**
632         * @param metaDataFilter the metaDataFilter to set
633         */
634        public void setMetaDataFilter(LinkedHashMap<String, String> metaDataFilter) {
635                this.metaDataFilter = metaDataFilter;
636        }
637
638        /**
639         * @return the coefficientsList
640         */
641        public LinkedHashMap<String, CoxCoefficient> getCoefficientsList() {
642                return coefficientsList;
643        }
644
645        /**
646         * @return the waldTestInfo
647         */
648        public WaldTestInfo getWaldTestInfo() {
649                return waldTestInfo;
650        }
651
652        /**
653         * @return the imat
654         */
655        public double[][] getImat() {
656                return imat;
657        }
658
659        /**
660         * @return the naive_imat
661         */
662        public double[][] getNaive_imat() {
663                return naive_imat;
664        }
665
666        /**
667         * @param degreeFreedom the degreeFreedom to set
668         */
669        public void setDegreeFreedom(double degreeFreedom) {
670                this.degreeFreedom = degreeFreedom;
671        }
672
673        /**
674         * @param waldTestInfo the waldTestInfo to set
675         */
676        public void setWaldTestInfo(WaldTestInfo waldTestInfo) {
677                this.waldTestInfo = waldTestInfo;
678        }
679}