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}