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 java.util.ArrayList;
024
025/**
026 *
027 * @author Scooter Willis <willishf at gmail dot com>
028 */
029public class CoxScore {
030
031        /**
032         *
033         * @param method
034         * @param survivalInfoList
035         * @param coxInfo
036         * @param useStrata
037         * @return
038         */
039        public static double[][] process(CoxMethod method, ArrayList<SurvivalInfo> survivalInfoList, CoxInfo coxInfo, boolean useStrata) {
040                int i, j, k;
041                double temp;
042                int n = survivalInfoList.size();
043
044                ArrayList<String> variables = new ArrayList<String>(coxInfo.getCoefficientsList().keySet());
045                int nvar = variables.size();
046
047                double deaths;
048                int dd;
049                double[] time = new double[n];
050                double[] status = new double[n];
051                double[] strata = new double[n];
052                double[] weights = new double[n];
053                double[] score = new double[n];
054                double[] a = new double[nvar];
055                double[] a2 = new double[nvar];
056                double denom = 0, e_denom;
057                double risk;
058                double[][] covar = new double[nvar][n];
059                double[][] resid = new double[nvar][n];
060                double hazard, meanwt;
061                double downwt, temp2;
062                double mean;
063
064                //  n = *nx;
065                //  nvar  = *nvarx;
066                for (int p = 0; p < n; p++) {
067                        SurvivalInfo si = survivalInfoList.get(p);
068                        time[p] = si.getTime();
069                        status[p] = si.getStatus();
070                        if (useStrata) {
071                                strata[p] = si.getStrata();
072                        } else {
073                                strata[p] = 0;
074                        }
075                        weights[p] = si.getWeight();
076                        score[p] = si.getScore();
077
078                        for(int v = 0; v < variables.size(); v++){
079                                String variable = variables.get(v);
080                                Double value = si.getVariable(variable);
081                                covar[v][p] = value;
082                        }
083
084                }
085
086
087
088                //  a = scratch;
089                //  a2 = a+nvar;
090        /*
091                 **  Set up the ragged array
092                 */
093                //   covar=  dmatrix(covar2, n, nvar);
094                //   resid=  dmatrix(resid2, n, nvar);
095
096                e_denom = 0;
097                deaths = 0;
098                meanwt = 0;
099                for (i = 0; i < nvar; i++) {
100                        a2[i] = 0;
101                }
102                strata[n - 1] = 1;  /*failsafe */
103                for (i = n - 1; i >= 0; i--) {
104                        if (strata[i] == 1) {
105                                denom = 0;
106                                for (j = 0; j < nvar; j++) {
107                                        a[j] = 0;
108                                }
109                        }
110
111                        risk = score[i] * weights[i];
112                        denom += risk;
113                        if (status[i] == 1) {
114                                deaths++;
115                                e_denom += risk;
116                                meanwt += weights[i];
117                                for (j = 0; j < nvar; j++) {
118                                        a2[j] += risk * covar[j][i];
119                                }
120                        }
121                        for (j = 0; j < nvar; j++) {
122                                a[j] += risk * covar[j][i];
123                                resid[j][i] = 0;
124                        }
125
126                        if (deaths > 0 && (i == 0 || strata[i - 1] == 1 || time[i] != time[i - 1])) {
127                                /* last obs of a set of tied death times */
128                                if (deaths < 2 || method == CoxMethod.Breslow) {
129                                        hazard = meanwt / denom;
130                                        for (j = 0; j < nvar; j++) {
131                                                temp = (a[j] / denom);     /* xbar */
132                                                for (k = i; k < n; k++) {
133                                                        temp2 = covar[j][k] - temp;
134                                                        if (time[k] == time[i] && status[k] == 1) {
135                                                                resid[j][k] += temp2;
136                                                        }
137                                                        resid[j][k] -= temp2 * score[k] * hazard;
138                                                        if (strata[k] == 1) {
139                                                                break;
140                                                        }
141                                                }
142                                        }
143                                } else {  /* the harder case */
144                                        meanwt /= deaths;
145                                        for (dd = 0; dd < deaths; dd++) {
146                                                downwt = dd / deaths;
147                                                temp = denom - downwt * e_denom;
148                                                hazard = meanwt / temp;
149                                                for (j = 0; j < nvar; j++) {
150                                                        mean = (a[j] - downwt * a2[j]) / temp;
151                                                        for (k = i; k < n; k++) {
152                                                                temp2 = covar[j][k] - mean;
153                                                                if (time[k] == time[i] && status[k] == 1) {
154                                                                        resid[j][k] += temp2 / deaths;
155                                                                        resid[j][k] -= temp2 * score[k] * hazard
156                                                                                        * (1 - downwt);
157                                                                } else {
158                                                                        resid[j][k] -= temp2 * score[k] * hazard;
159                                                                }
160                                                                if (strata[k] == 1) {
161                                                                        break;
162                                                                }
163                                                        }
164                                                }
165                                        }
166                                }
167                                e_denom = 0;
168                                deaths = 0;
169                                meanwt = 0;
170                                for (j = 0; j < nvar; j++) {
171                                        a2[j] = 0;
172                                }
173                        }
174                }
175
176                for (int p = 0; p < n; p++) {
177                        SurvivalInfo si = survivalInfoList.get(p);
178                        for (int v = 0; v < variables.size(); v++) {
179                                si.setResidualVariable(variables.get(v), resid[v][p]);
180                        }
181
182                }
183
184                //appears to be backward internally
185                double[][] flipresid = new double[n][nvar];
186
187                for(int s = 0; s < resid.length; s++){
188                        for(int t = 0; t  < resid[0].length; t++){
189                                flipresid[t][s] = resid[s][t];
190                        }
191                }
192
193                return flipresid;
194
195        }
196
197        /**
198         * @param args the command line arguments
199         */
200        public static void main(String[] args) {
201                // TODO code application logic here
202        }
203}