org.biojava.nbio.survival.cox.ResidualsCoxph Maven / Gradle / Ivy
The newest version!
/*
* BioJava development code
*
* This code may be freely distributed and modified under the
* terms of the GNU Lesser General Public Licence. This should
* be distributed with the code. If you do not have a copy,
* see:
*
* http://www.gnu.org/copyleft/lesser.html
*
* Copyright for this code is held jointly by the individual
* authors. These should be listed in @author doc comments.
*
* For more information on the BioJava project and its aims,
* or to join the biojava-l mailing list, visit the home page
* at:
*
* http://www.biojava.org/
*
*/
package org.biojava.nbio.survival.cox;
import org.biojava.nbio.survival.cox.matrix.Matrix;
import java.util.ArrayList;
import java.util.LinkedHashMap;
//import java.util.Collections;
/**
*
* @author Scooter Willis
*/
public class ResidualsCoxph {
/**
*
*/
public enum Type {
/**
*
*/
dfbeta,
/**
*
*/
dfbetas,
/**
*
*/
score;
}
/**
*
* @param ci
* @param type
* @param useWeighted
* @param cluster
* @return
* @throws Exception
*/
public static double[][] process(CoxInfo ci, Type type, boolean useWeighted, ArrayList cluster) throws Exception {
Type otype = type;
if (type == Type.dfbeta || type == Type.dfbetas) {
type = Type.score;
//if missing weighted is a required so never missing
} //64 2 625 310
double[][] rr = null;
if (type == Type.score) {
rr = CoxScore.process(ci.method, ci.survivalInfoList, ci, false);
}
//debug
// if (false) {
// for (int i = 0; i < ci.survivalInfoList.size(); i++) {
// SurvivalInfo si = ci.survivalInfoList.get(i);
// System.out.print("residuals " + si.getOrder() + " " + si.getClusterValue());
// for (int j = 0; j < 2; j++) {
// System.out.print(" " + rr[i][j]);
// }
// System.out.println();
// }
// }
double[][] vv = null;
if (ci.getNaiveVariance() != null) {
vv = ci.getNaiveVariance();
} else {
vv = ci.getVariance();
}
if (otype == Type.dfbeta) {
//rr <- rr %*% vv
rr = Matrix.multiply(rr, vv);
} else if (otype == Type.dfbetas) {
//rr <- (rr %*% vv) %*% diag(sqrt(1/diag(vv)))
double[][] d1 = Matrix.multiply(rr, vv);
double[][] d2 = Matrix.diag(Matrix.sqrt(Matrix.oneDivide(Matrix.diag(vv))));
rr = Matrix.multiply(d1, d2);
}
if (useWeighted) {
double[] weighted = ci.getWeighted();
rr = Matrix.scale(rr, weighted);
}
if (cluster != null && cluster.size() > 0) {
rr = rowsum(rr, cluster);
}
return rr;
}
/**
* From R in residuals.coxph.S rowsum(rr, collapse)
*
* @param rr
* @param sets
* @return
*/
private static double[][] rowsum(double[][] rr, ArrayList sets) throws Exception {
LinkedHashMap sumMap = new LinkedHashMap<>();
if (rr.length != sets.size()) {
throw new Exception("Cluster value for each sample are not of equal length n=" + rr.length + " cluster length=" + sets.size());
}
double[][] sum = null;
for (int j = 0; j < rr[0].length; j++) {
for (int i = 0; i < sets.size(); i++) {
String s = sets.get(i);
Double v = sumMap.get(s); //get in order
if (v == null) {
v = 0.0;
}
v = v + rr[i][j];
sumMap.put(s, v);
}
if (sum == null) {
sum = new double[sumMap.size()][rr[0].length];
}
ArrayList index = new ArrayList<>(sumMap.keySet());
//sorting does seem to make a difference in test cases at the .0000000001
// ArrayList in = new ArrayList();
// for (String s : index) {
// in.add(Integer.parseInt(s));
// }
// Collections.sort(index);
for (int m = 0; m < index.size(); m++) {
String key = index.get(m);
sum[m][j] = sumMap.get(key);
}
sumMap.clear();
}
return sum;
}
/**
* @param args the command line arguments
*/
public static void main(String[] args) {
// TODO code application logic here
}
}