[Genabel-commits] r963 - in branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3: . src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 30 16:52:56 CEST 2012
Author: lckarssen
Date: 2012-09-30 16:52:56 +0200 (Sun, 30 Sep 2012)
New Revision: 963
Added:
branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/coxscore.c
Modified:
branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/Makefile
branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/data.h
branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/main.cpp
branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/reg1.h
Log:
In ProbABEL-pacox-robust branch:
Applied patch by Alisa K. Manning.
Work in progress.
Modified: branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/Makefile
===================================================================
--- branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/Makefile 2012-09-30 14:37:12 UTC (rev 962)
+++ branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/Makefile 2012-09-30 14:52:56 UTC (rev 963)
@@ -11,7 +11,7 @@
## in this tree only build pacoxph. Leave these out: $(LOGREG) $(LINREG)
REGFILES = $(SRCDIR)/data.h $(SRCDIR)/mematrix.h $(SRCDIR)/reg1.h $(SRCDIR)/usage.h $(SRCDIR)/main.cpp
-COXBASE = $(SRCDIR)/chinv2 $(SRCDIR)/cholesky2 $(SRCDIR)/chsolve2 $(SRCDIR)/dmatrix $(SRCDIR)/coxfit2
+COXBASE = $(SRCDIR)/chinv2 $(SRCDIR)/cholesky2 $(SRCDIR)/chsolve2 $(SRCDIR)/dmatrix $(SRCDIR)/coxfit2 $(SRCDIR)/coxscore
COXSRC = $(COXBASE:=.c)
COXOBJ = $(COXBASE:=.o)
Added: branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/coxscore.c
===================================================================
--- branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/coxscore.c (rev 0)
+++ branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/coxscore.c 2012-09-30 14:52:56 UTC (rev 963)
@@ -0,0 +1,127 @@
+/* $Id: coxscore.c 11357 2009-09-04 15:22:46Z therneau $
+**
+** Compute the score residuals for a Cox model
+**
+** Input
+** nx number of subjects
+** nvarx number of variables in the covariance matrix
+** y matrix of time and status values
+** strata =1 for the last obs of each strata
+** covar2 the matrix of covariates, rows=variables, columns=subjects
+** (the S executive stores matrices in the Fortran ordering)
+** score the vector of subject scores, i.e., exp(beta*z)
+** weights case weight
+** method ==1 for efron method
+**
+** Output
+** resid a matrix of the same shape as x
+**
+** Scratch
+** scratch, from which a and a2 are carved
+**
+** Data must be sorted by strata, ascending time within strata, death before
+** censor within time.
+*/
+#include <stdio.h>
+#include "survS.h"
+#include "survproto.h"
+
+void coxscore(Sint *nx, Sint *nvarx, double *y,
+ double *covar2, Sint *strata, double *score,
+ double *weights, Sint *method, double *resid2,
+ double *scratch)
+ {
+ int i,j, k;
+ double temp;
+ int n, nvar;
+ double deaths;
+ int dd;
+ double *time, *status;
+ double *a, *a2;
+ double denom=0, e_denom;
+ double risk;
+ double **covar;
+ double **resid;
+ double hazard, meanwt;
+ double downwt, temp2;
+ double mean;
+
+ n = *nx;
+ nvar = *nvarx;
+ time = y;
+ status = y+n;
+ a = scratch;
+ a2 = a+nvar;
+ /*
+ ** Set up the ragged array
+ */
+ covar= dmatrix(covar2, n, nvar);
+ resid= dmatrix(resid2, n, nvar);
+
+ e_denom=0;
+ deaths=0;
+ meanwt=0;
+ for (i=0; i<nvar; i++) a2[i] =0;
+ strata[n-1] =1; /*failsafe */
+ for (i=n-1; i >=0; i--) {
+ if (strata[i]==1) {
+ denom =0;
+ for (j=0; j<nvar; j++) a[j] =0;
+ }
+
+ risk = score[i] * weights[i];
+ denom += risk;
+ if (status[i]==1) {
+ deaths++;
+ e_denom += risk;
+ meanwt += weights[i];
+ for (j=0; j<nvar; j++) a2[j] += risk*covar[j][i];
+ }
+ for (j=0; j<nvar; j++) {
+ a[j] += risk * covar[j][i];
+ resid[j][i] =0;
+ }
+
+ if (deaths>0 && (i==0 || strata[i-1]==1 || time[i]!=time[i-1])){
+ /* last obs of a set of tied death times */
+ if (deaths <2 || *method==0) {
+ hazard = meanwt/denom;
+ for (j=0; j<nvar; j++) {
+ temp = (a[j]/denom); /* xbar */
+ for (k=i; k<n; k++) {
+ temp2 = covar[j][k] - temp;
+ if (time[k]==time[i] && status[k]==1)
+ resid[j][k] += temp2;
+ resid[j][k] -= temp2* score[k] * hazard;
+ if (strata[k]==1) break;
+ }
+ }
+ }
+ else { /* the harder case */
+ meanwt /= deaths;
+ for (dd=0; dd<deaths; dd++) {
+ downwt = dd/deaths;
+ temp = denom - downwt* e_denom;
+ hazard = meanwt/temp;
+ for (j=0; j<nvar; j++) {
+ mean = (a[j] - downwt*a2[j])/ temp;
+ for (k=i; k<n; k++) {
+ temp2 = covar[j][k] - mean;
+ if (time[k]==time[i] && status[k]==1) {
+ resid[j][k] += temp2/deaths;
+ resid[j][k] -= temp2 * score[k] * hazard *
+ (1 - downwt);
+ }
+ else resid[j][k]-= temp2*score[k] * hazard;
+ if (strata[k]==1) break;
+ }
+ }
+ }
+ }
+ e_denom =0;
+ deaths =0;
+ meanwt =0;
+ for (j=0; j<nvar; j++) a2[j] =0;
+ }
+ }
+ }
Modified: branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/data.h
===================================================================
--- branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/data.h 2012-09-30 14:37:12 UTC (rev 962)
+++ branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/data.h 2012-09-30 14:52:56 UTC (rev 963)
@@ -405,6 +405,7 @@
mematrix<int> strata;
mematrix<double> X;
mematrix<int> order;
+ mematrix<double> Y; //AKM - call to coxscore function requires response to be stored in a Nx2 matrix
coxph_data(phedata &phed, gendata &gend, int snpnum)
{
@@ -427,6 +428,7 @@
offset.reinit(nids,1);
strata.reinit(nids,1);
order.reinit(nids,1);
+ Y.reinit(nids,2); // AKM
for (int i=0;i<nids;i++)
{
// X.put(1.,i,0);
@@ -437,6 +439,8 @@
fprintf(stderr,"coxph_data: status not 0/1 (right order: id, fuptime, status ...)\n",phed.noutcomes);
exit(1);
}
+ Y.put((phed.Y).get(i,0),i,0); // event time
+ Y.put((phed.Y).get(i,1),i,1); // event 1=yes 0=no
}
for (int j=0;j<phed.ncov;j++)
for (int i=0;i<nids;i++)
@@ -483,6 +487,7 @@
offset = reorder(offset,order);
X = reorder(X,order);
X = transpose(X);
+ Y = reorder(Y,order); // AKM
// X.print();
// offset.print();
// weights.print();
Modified: branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/main.cpp
===================================================================
--- branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/main.cpp 2012-09-30 14:37:12 UTC (rev 962)
+++ branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/main.cpp 2012-09-30 14:52:56 UTC (rev 963)
@@ -294,7 +294,10 @@
#if COXPH
if(inverse_filename != NULL) {std::cerr<<"ERROR: mmscore is forbidden for cox regression\n";exit(1);}
- if (robust) {std::cerr<<"ERROR: robust standard errors not implemented for Cox regression (drop us e-mail if you really need that)\n";exit(1);}
+ if (robust) {
+ //std::cerr<<"ERROR: robust standard errors not implemented for Cox regression (drop us e-mail if you really need that)\n";exit(1);
+ cout << "Robust standard errors implemented for Cox regression\n";
+ }
#endif
@@ -337,7 +340,7 @@
nrd.estimate(nrgd,0,CHOLTOL,0, interaction, ngpreds, invvarmatrix, robust, 1);
#elif COXPH
coxph_reg nrd(nrgd);
- nrd.estimate(nrgd,0,MAXITER,EPS,CHOLTOL,0, interaction, ngpreds, 1);
+ nrd.estimate(nrgd,0,MAXITER,EPS,CHOLTOL,0, interaction, ngpreds, 1, robust);
#endif
null_loglik = nrd.loglik;
@@ -427,10 +430,10 @@
//Han Chen
*outfile[0] << sep << "beta_SNP_A1A2_" << phd.model_terms[interaction_cox] << sep << "sebeta_SNP_A1A2_" << phd.model_terms[interaction_cox]
<< sep << "beta_SNP_A1A1_" << phd.model_terms[interaction_cox] << sep << "sebeta_SNP_A1A1_" << phd.model_terms[interaction_cox];
- #if !COXPH
+ // #if !COXPH
if(inverse_filename == NULL && !allcov) *outfile[0] << sep << "cov_SNP_A1A2_int_SNP_" << phd.model_terms[interaction_cox]
<< sep << "cov_SNP_A1A1_int_SNP_" << phd.model_terms[interaction_cox];
- #endif
+ // #endif
//Oct 26, 2009
for (int file=1; file<outfile.size() ; file++)
{
@@ -496,9 +499,9 @@
if(inverse_filename == NULL)
//Han Chen
{
- #if !COXPH
+ // #if !COXPH
if(interaction != 0 && !allcov) *outfile[0] << sep << "cov_SNP_int_SNP_" << phd.model_terms[interaction_cox];
- #endif
+ // #endif
*outfile[0] << sep << "chi2_SNP";
}
//Oct 26, 2009
@@ -672,7 +675,7 @@
}
#elif COXPH
coxph_reg rd(rgd);
- rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model, interaction, true, ngpreds);
+ rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model, interaction, true, ngpreds,0,robust);
#endif
if(!allcov && model==0 && interaction==0) start_pos=rd.beta.nrow-2;
@@ -686,7 +689,7 @@
{
*beta_sebeta[model] << sep << rd.beta[pos] << sep << rd.sebeta[pos];
//Han Chen
- #if !COXPH
+// #if !COXPH
if (inverse_filename == NULL && !allcov && interaction != 0)
{
if (pos>start_pos)
@@ -698,7 +701,7 @@
{*covvalue[model] << rd.covariance[pos-1];}
}
}
- #endif
+// #endif
//Oct 26, 2009
}
@@ -734,14 +737,14 @@
*beta_sebeta[model] << sep << "nan" << sep << "nan";
}
//Han Chen
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{if (model==0)
{*covvalue[model] << "nan" << sep << "nan";}
else
{*covvalue[model] << "nan";}
}
- #endif
+// #endif
//Oct 26, 2009
*chi2[model] << "nan";
}
@@ -752,44 +755,44 @@
//Han Chen
*outfile[0] << beta_sebeta[0]->str() << sep;
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{
*outfile[0] << covvalue[0]->str() << sep;
}
- #endif
+// #endif
*outfile[0] << chi2[0]->str() << "\n";
*outfile[1] << beta_sebeta[1]->str() << sep;
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{
*outfile[1] << covvalue[1]->str() << sep;
}
- #endif
+// #endif
*outfile[1] << chi2[1]->str() << "\n";
*outfile[2] << beta_sebeta[2]->str() << sep;
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{
*outfile[2] << covvalue[2]->str() << sep;
}
- #endif
+// #endif
*outfile[2] << chi2[2]->str() << "\n";
*outfile[3] << beta_sebeta[3]->str() << sep;
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{
*outfile[3] << covvalue[3]->str() << sep;
}
- #endif
+// #endif
*outfile[3] << chi2[3]->str() << "\n";
*outfile[4] << beta_sebeta[4]->str() << sep;
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{
*outfile[4] << covvalue[4]->str() << sep;
}
- #endif
+// #endif
*outfile[4] << chi2[4]->str() << "\n";
//Oct 26, 2009
@@ -832,7 +835,7 @@
}
#elif COXPH
coxph_reg rd(rgd);
- rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model, interaction, true, ngpreds);
+ rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model, interaction, true, ngpreds,0,robust);
#endif
if(!allcov && interaction==0) start_pos=rd.beta.nrow-1;
@@ -845,12 +848,12 @@
{
*beta_sebeta[0] << sep << rd.beta[pos] << sep << rd.sebeta[pos];
//Han Chen
- #if !COXPH
+// #if !COXPH
if (inverse_filename == NULL && !allcov && interaction != 0)
{if (pos>start_pos)
{*covvalue[0] << rd.covariance[pos-1];}
}
- #endif
+// #endif
//Oct 26, 2009
}
@@ -893,10 +896,10 @@
if(inverse_filename == NULL)
{
//Han Chen
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{*covvalue[0] << "nan";}
- #endif
+// #endif
//Oct 26, 2009
*chi2[0] << "nan";
}
@@ -906,10 +909,10 @@
{
//Han Chen
*outfile[0] << beta_sebeta[0]->str() << sep;
- #if !COXPH
+// #if !COXPH
if (!allcov && interaction !=0)
{*outfile[0] << covvalue[0]->str() << sep;}
- #endif
+// #endif
*outfile[0] << chi2[model]->str() << "\n";
//Oct 26, 2009
}
Modified: branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/reg1.h
===================================================================
--- branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/reg1.h 2012-09-30 14:37:12 UTC (rev 962)
+++ branches/ProbABEL-pacox-robust/ProbABEL-pacox/v.0.1-3/src/reg1.h 2012-09-30 14:52:56 UTC (rev 963)
@@ -785,6 +785,12 @@
double *work, double *eps, double *tol_chol,
double *sctest);
+void coxscore(int *nx, int *nvarx, double *y,
+ double *covar2, int *strata, double *score,
+ double *weights, int *method, double *resid2,
+ double *scratch);
+
+
class coxph_reg
{
public:
@@ -795,7 +801,10 @@
double loglik;
double chi2_score;
int niter;
+
+ mematrix<double> covariance; // AKM
+
coxph_reg(coxph_data &cdata)
{
beta.reinit(cdata.X.nrow,1);
@@ -804,13 +813,22 @@
sigma2=-1.;
chi2_score=-1.;
niter=0;
+
+ // AKM: This array stores the covariance between the last term in the
+ // regression model and all previous terms. Used for interaction analysis to
+ // store the covariance(s) between SNPxE and SNP OR ((SNPAAxE and SNPAA) and (SNPABxE and SNPAB))
+ covariance.reinit(cdata.X.nrow-1,1);
+
+ residuals.reinit(cdata.X.ncol,1); // AKM - residuals
+
+
}
~coxph_reg()
{
// delete beta;
// delete sebeta;
}
- void estimate(coxph_data &cdata, int verbose, int maxiter, double eps, double tol_chol, int model, int interaction, int ngpreds, bool iscox, int nullmodel=0)
+ void estimate(coxph_data &cdata, int verbose, int maxiter, double eps, double tol_chol, int model, int interaction, int ngpreds, bool iscox, int nullmodel=0,int robust=0)
{
// cout << "model = " << model << "\n";
// cdata.X.print();
@@ -819,6 +837,15 @@
int length_beta = X.nrow;
beta.reinit(length_beta,1);
sebeta.reinit(length_beta,1);
+ // AKM
+ if (length_beta>1) {
+ if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2) {
+ covariance.reinit(length_beta-2,1);
+ } else {
+ covariance.reinit(length_beta-1,1);
+ }
+ }
+
mematrix<double> newoffset = cdata.offset;
newoffset = cdata.offset - (cdata.offset).column_mean(0);
mematrix<double> means(X.nrow,1);
@@ -837,7 +864,94 @@
imat.data,loglik_int,&flag,
work,&eps,&tol_chol,
&sctest);
- for (int i=0;i<X.nrow;i++) sebeta[i]=sqrt(imat.get(i,i));
+
+ // AKM - robust variance calculation
+ if(robust == 1) {
+
+ mematrix<double> Xorig = t_apply_model(cdata.X,model, interaction, ngpreds, iscox, nullmodel);
+ mematrix<double> score(X.ncol,1);
+
+ mematrix<double> tX = transpose(X);
+
+ double zbeta, risk;
+
+ for (int person=0; person<cdata.nids; person++) {
+ // means.data is the offset of the linear predictor
+ zbeta = 0; /* form the term beta*z (vector mult) */
+
+ for (int i=0; i<tX.ncol; i++) {
+ zbeta += beta[i]*tX.get(person,i);
+
+ }
+ double temp = exp(zbeta);
+
+ score[person] = temp;
+ }
+
+
+ // call coxscore function (some of the same inputs as coxfit2 above
+ // Returned is a matrix the same dimension as X (need temp matrix!)
+ //mematrix<double> tempresiduals(X.nrow,X.ncol);
+ //mematrix<double> tempa(2,tX.ncol);
+ int method = 1;
+
+ mematrix<double> tempweights(1,cdata.nids);
+ tempweights.reinit(cdata.nids,1);
+ //tempresiduals.reinit(X.nrow,X.ncol);
+ double tempresiduals[X.nrow*X.ncol];
+ double tempa[2*tX.ncol];
+
+ for(int i=0;i<cdata.nids;i++){
+ tempweights[i] = 1;
+ }
+
+ coxscore(&cdata.nids,&X.nrow,cdata.Y.data,Xorig.data,
+ cdata.strata.data,score.data,tempweights.data,&method,
+ tempresiduals, // data is returned here.
+ tempa); // a is a matrix two rows and nvar columns
+ cout << "PRINTING RESULT\n";
+ for(int k=0;k<X.nrow*X.ncol;k++) {
+ cout << "\t" << tempresiduals[k];
+ }
+ cout << "\nDONE PRINTING RESULT\n";
+
+// need to put values of tempresiduals into a mematrix -- TO DO
+// residuals = transpose(tempresiduals)*imat;
+
+ // must multiply by weight?? weight is a matrix of all 1's so unnecessary
+
+ if(1) {
+ cout << "RESIDUALS\n";
+// residuals.print();
+ }
+
+// imat = transpose(residuals)*residuals;
+ //cdata.residuals.print();
+
+ cout << "END ROBUST\n\n";
+ }
+
+
+ imat.print();
+
+ for (int i=0;i<X.nrow;i++) {
+ sebeta[i]=sqrt(imat.get(i,i));
+ if (i>0) { // not first variable
+ if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2) {
+ // SNPAA and SNPAB in model
+ if (i>1) {
+ double covval=imat.get(i,i-2); // grab cov BETWEEN SNPAA*E and SNPAA, for example
+ covariance.put(covval,i-2,0);
+ }
+ } else { // one term in model for SNP
+
+ double covval=imat.get(i,i-1);
+ covariance.put(covval,i-1,0);
+ }
+ }
+ }
+
+// for (int i=0;i<X.nrow;i++) sebeta[i]=sqrt(imat.get(i,i));
loglik = loglik_int[1];
niter = maxiter;
}
More information about the Genabel-commits
mailing list