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 */
021/*
022 * /* Ported from $Id: coxmart.c 11166 2008-11-24 22:10:34Z therneau $
023 *
024 ** Compute the martingale residual for a Cox model
025 **
026 ** Input
027 **      n       number of subjects
028 **      method  will be ==1 for the Efron method
029 **      time    vector of times
030 **      status  vector of status values
031 **      score   the vector of subject scores, i.e., exp(beta*z)
032 **      strata  is =1 for the last obs of a strata
033 **      mark    carried forward from the coxfit routine
034 **
035 ** Output
036 **      expected the expected number of events for the subject
037 **
038 ** The martingale residual is more of a nuisance for the Efron method
039 **
040 */
041package org.biojava.nbio.survival.cox;
042
043import java.util.ArrayList;
044
045/**
046 *
047 * @author Scooter Willis <willishf at gmail dot com>
048 */
049public class CoxMart {
050
051        /**
052         *
053         * @param method
054         * @param survivalInfoList
055         * @param useStrata
056         * @return
057         */
058        static public double[] process(CoxMethod method, ArrayList<SurvivalInfo> survivalInfoList, boolean useStrata) {
059                int i, j;
060                int lastone;
061                int n = survivalInfoList.size();
062                double deaths, denom = 0, e_denom = 0;
063                double hazard;
064                double temp, wtsum;
065                double downwt;
066                double[] time = new double[n];
067                double[] status = new double[n];
068                double[] strata = new double[n];
069                double[] wt = new double[n];
070                double[] score = new double[n];
071
072                double[] expect = new double[survivalInfoList.size()];
073
074                for (int p = 0; p < n; p++) {
075                        SurvivalInfo si = survivalInfoList.get(p);
076                        time[p] = si.getTime();
077                        status[p] = si.getStatus();
078                        if (useStrata) {
079                                strata[p] = si.getStrata();
080                        } else {
081                                strata[p] = 0;
082                        }
083                        wt[p] = si.getWeight();
084                        score[p] = si.getScore();
085                }
086
087                strata[n - 1] = 1;  /*failsafe */
088                /* Pass 1-- store the risk denominator in 'expect' */
089                for (i = n - 1; i >= 0; i--) {  // Error because of no bounds checking in C it is an error on the get i - 1
090                 //   SurvivalInfo si = survivalInfoList.get(i);
091
092                        if (strata[i] == 1) {
093                                denom = 0; //strata[i]
094                        }
095                        denom += score[i] * wt[i];      //score[i]*wt[i];
096                        if (i == 0 || strata[i - 1] == 1 || time[i - 1] != time[i]) //strata[i-1]==1 ||  time[i-1]!=time[i]
097                        {
098                         //   si.setResidual(denom);
099                                expect[i] = denom;
100                        } else {
101                         //   si.setResidual(0); //expect[i] =0;
102                                expect[i] = 0;
103                        }
104                }
105
106                /* Pass 2-- now do the work */
107                deaths = 0;
108                wtsum = 0;
109                e_denom = 0;
110                hazard = 0;
111                lastone = 0;
112                for (i = 0; i < n; i++) {
113        //         SurvivalInfo si = survivalInfoList.get(i);
114        //         SurvivalInfo sip1 = null;
115        //         if (i + 1 < n) {
116        //             sip1 = survivalInfoList.get(i + 1);
117        //         }
118        //         if (si.getResidual() != 0) {
119        //             denom = si.getResidual();
120        //         }
121                        if(expect[i] != 0)
122                                denom = expect[i];
123        //         si.setResidual(status[i]);//expect[i] = status[i];
124                        expect[i] = status[i];
125
126                        deaths += status[i]; //status[i];
127                        wtsum += status[i] * wt[i]; // status[i]*wt[i];
128                        e_denom += score[i] * status[i] * wt[i];//score[i]*status[i] *wt[i];
129                        if ( strata[i] == 1 || time[i + 1] != time[i]) { //strata[i]==1 ||  time[i+1]!=time[i]
130                /*last subject of a set of tied times */
131                                if (deaths < 2 || method == CoxMethod.Breslow) { //*method==0
132                                        hazard += wtsum / denom;
133                                        for (j = lastone; j <= i; j++) {
134                                        //     SurvivalInfo sj = survivalInfoList.get(j);
135                                        //     double res =  sj.getResidual() - score[j] * hazard;
136                                                expect[j] -= score[j] * hazard;
137                                        //     sj.setResidual(res); //expect[j] -= score[j]*hazard;
138
139                                        }
140                                } else {
141                                        temp = hazard;
142                                        wtsum /= deaths;
143                                        for (j = 0; j < deaths; j++) {
144                                                downwt = j / deaths;
145                                                hazard += wtsum / (denom - e_denom * downwt);
146                                                temp += wtsum * (1 - downwt) / (denom - e_denom * downwt);
147                                        }
148                                        for (j = lastone; j <= i; j++) {
149                                         //   SurvivalInfo sj = survivalInfoList.get(j);
150                                                if (status[j] == 0) {
151                                                        //this appears to be an error for = - versus -=
152                                                 //   double res = -score[j] * hazard;
153                                                        expect[j] = -score[j] * hazard;
154                                                 //   sj.setResidual(res);//expect[j] = -score[j]*hazard; This appears to be an error of -score vs -=
155                                                } else {
156                                                        // double res = sj.getResidual() - score[j] * temp;
157                                                        expect[j] -= score[j] * temp;
158                                                        //expect[j] -=  score[j]* temp;
159                                                }
160                                        }
161                                }
162                                lastone = i + 1;
163                                deaths = 0;
164                                wtsum = 0;
165                                e_denom = 0;
166                        }
167                        if (strata[i] == 1) {
168                                hazard = 0;
169                        }
170                }
171
172
173                for (j = lastone; j < n; j++) {
174                        expect[j] -= score[j] * hazard;
175                }
176
177                return expect;
178        }
179
180        /**
181         * @param args the command line arguments
182         */
183        public static void main(String[] args) {
184                // TODO code application logic here
185        }
186}