[Genabel-commits] r916 - in tags/ProbABEL: . v.0.1-3/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 11 22:48:47 CEST 2012
Author: lckarssen
Date: 2012-06-11 22:48:47 +0200 (Mon, 11 Jun 2012)
New Revision: 916
Added:
tags/ProbABEL/v.0.2.0/
Modified:
tags/ProbABEL/v.0.1-3/src/data.h
tags/ProbABEL/v.0.1-3/src/reg1.h
Log:
Prepare for tagging ProbABEL v0.2.0 release
Modified: tags/ProbABEL/v.0.1-3/src/data.h
===================================================================
--- tags/ProbABEL/v.0.1-3/src/data.h 2012-06-11 07:45:12 UTC (rev 915)
+++ tags/ProbABEL/v.0.1-3/src/data.h 2012-06-11 20:48:47 UTC (rev 916)
@@ -19,19 +19,19 @@
}
char tmp[100];
- for (int i=0;i<nphenocols;i++)
+ for (int i=0;i<nphenocols;i++)
{
fscanf(infile,"%s",&tmp);
// printf("%s ",tmp);
} //printf("\n");
- unsigned short int * allmeasured = new unsigned short int [npeople];
+ unsigned short int * allmeasured = new unsigned short int [npeople];
int nids = 0;
for (int i = 0;i<npeople;i++)
{
allmeasured[i] = 1;
fscanf(infile,"%s",&tmp);
- for (int j=1;j<nphenocols;j++)
+ for (int j=1;j<nphenocols;j++)
{
fscanf(infile,"%s",&tmp);
if (tmp[0]=='N' || tmp[0]=='n') allmeasured[i]=0;
@@ -42,7 +42,7 @@
return(nids);
}
-class phedata
+class phedata
{
public:
int nids_all;
@@ -56,7 +56,7 @@
std::string model;
std::string * model_terms;
int n_model_terms;
- phedata(char * fname, int noutc, int npeople, int interaction, bool iscox)
+ phedata(char * fname, int noutc, int npeople, int interaction, bool iscox)
{
static const unsigned int BFS = 1000;
std::ifstream myfile(fname);
@@ -93,18 +93,18 @@
};
myfile.close();
}
- else
+ else
{
- fprintf(stderr,"Unable to open file %s\n",fname);
+ fprintf(stderr,"Unable to open file %s\n",fname);
exit(1);
}
fprintf(stdout,"Actual number of people in phenofile = %d",npeople);
- if (savenpeople>0)
+ if (savenpeople>0)
{
npeople = savenpeople;
fprintf(stdout,"; using only %d first\n",npeople);
}
- else
+ else
{
fprintf(stdout,"; using all of these\n");
}
@@ -124,7 +124,7 @@
model = "( ";
fscanf(infile,"%s",&tmp);
model = model + tmp;
- for (int i = 1;i < noutcomes;i++)
+ for (int i = 1;i < noutcomes;i++)
{
fscanf(infile,"%s",&tmp);
model = model + " , ";
@@ -138,15 +138,15 @@
model_terms[n_model_terms++] = "mu";
#endif
- if (nphenocols>noutcomes+1)
+ if (nphenocols>noutcomes+1)
{
fscanf(infile,"%s",&tmp);
model = model + tmp;
model_terms[n_model_terms++] = tmp;
- for (int i=(2+noutcomes);i<nphenocols;i++)
+ for (int i=(2+noutcomes);i<nphenocols;i++)
{
fscanf(infile,"%s",&tmp);
-
+
// if(iscox && ) {if(n_model_terms+1 == interaction-1) {continue;} }
// else {if(n_model_terms+1 == interaction) {continue;} }
model = model + " + ";
@@ -180,12 +180,12 @@
#endif
std::cout << "model: " << model << "\n";
- allmeasured = new unsigned short int [npeople];
+ allmeasured = new unsigned short int [npeople];
nids = 0;
for (int i = 0;i<npeople;i++)
{
allmeasured[i] = 1;
- for (int j=0;j<nphenocols;j++)
+ for (int j=0;j<nphenocols;j++)
{
fscanf(infile,"%s",&tmp);
if (j>0 && (tmp[0]=='N' || tmp[0]=='n')) allmeasured[i]=0;
@@ -201,14 +201,14 @@
idnames = new std::string [nids];
X.reinit(nids,ntmpcov);
Y.reinit(nids,noutcomes);
-
+
// second pass -- read the data
if ((infile=fopen(fname,"r"))==NULL) {
fprintf(stderr,"phedata: can not open file %s\n",fname);
exit(1);
}
- for (int i=0;i<nphenocols;i++)
+ for (int i=0;i<nphenocols;i++)
{
fscanf(infile,"%s",&tmp);
}
@@ -216,7 +216,7 @@
int k =0;
int m =0;
for (int i = 0;i<npeople;i++)
- if (allmeasured[i]==1)
+ if (allmeasured[i]==1)
{
fscanf(infile,"%s",&tmp);
idnames[m] = tmp;
@@ -225,14 +225,14 @@
fscanf(infile,"%s",&tmp);
Y.put(atof(tmp),m,j);
}
- for (int j=(1+noutcomes);j<nphenocols;j++)
+ for (int j=(1+noutcomes);j<nphenocols;j++)
{
fscanf(infile,"%s",&tmp);
X.put(atof(tmp),m,(j-1-noutcomes));
}
m++;
- }
- else
+ }
+ else
for (int j=0;j<nphenocols;j++) fscanf(infile,"%s",&tmp);
fclose(infile);
}
@@ -258,7 +258,7 @@
nsnps = insnps;
ngpreds = ingpreds;
int nids_all = npeople;
-
+
G.reinit(nids,(nsnps*ngpreds));
FILE * infile;
@@ -272,7 +272,7 @@
int k = 0;
for (int i = 0;i<npeople;i++)
- if (allmeasured[i]==1)
+ if (allmeasured[i]==1)
{
if (skipd>0)
{
@@ -294,10 +294,10 @@
exit(1);
}
}
- for (int j=1;j<skipd;j++) {
+ for (int j=1;j<skipd;j++) {
fscanf(infile,"%s",&tmp);
}
- for (int j=0;j<(nsnps*ngpreds);j++)
+ for (int j=0;j<(nsnps*ngpreds);j++)
{
int a = fscanf(infile,"%s",&tmp);
if (!a || a==EOF)
@@ -309,8 +309,8 @@
G.put(atof(tmp),k,j);
}
k++;
- }
- else
+ }
+ else
{
for (int j=0;j<skipd;j++) fscanf(infile,"%s",&tmp);
for (int j=0;j<(nsnps*ngpreds);j++) fscanf(infile,"%s",&tmp);
@@ -334,34 +334,34 @@
mematrix<double> X;
mematrix<double> Y;
- regdata(phedata &phed, gendata &gend, int snpnum)
+ regdata(phedata &phed, gendata &gend, int snpnum)
{
nids = gend.nids;
ngpreds = gend.ngpreds;
- if (snpnum>=0)
+ if (snpnum>=0)
ncov = phed.ncov + ngpreds;
- else
+ else
ncov = phed.ncov;
noutcomes = phed.noutcomes;
X.reinit(nids,(ncov+1));
Y.reinit(nids,noutcomes);
- for (int i=0;i<nids;i++)
+ for (int i=0;i<nids;i++)
{
X.put(1.,i,0);
Y.put((phed.Y).get(i,0),i,0);
}
- for (int j=1;j<=phed.ncov;j++)
- for (int i=0;i<nids;i++)
+ for (int j=1;j<=phed.ncov;j++)
+ for (int i=0;i<nids;i++)
X.put((phed.X).get(i,j-1),i,j);
- if (snpnum>0)
- for (int i=0;i<nids;i++)
- for (int j=0;j<ngpreds;j++)
+ if (snpnum>0)
+ for (int i=0;i<nids;i++)
+ for (int j=0;j<ngpreds;j++)
X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+1+j));
}
void update_snp(gendata &gend, int snpnum)
{
- for (int i=0;i<nids;i++)
- for (int j=0;j<ngpreds;j++)
+ for (int i=0;i<nids;i++)
+ for (int j=0;j<ngpreds;j++)
X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov+1-j-1));
// X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+1+(ngpreds-j+1)));
// X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+1+j));
@@ -395,109 +395,109 @@
class coxph_data
{
public:
- int nids;
- int ncov;
- int ngpreds;
- mematrix<double> weights;
- mematrix<double> stime;
- mematrix<int> sstat;
- mematrix<double> offset;
- mematrix<int> strata;
- mematrix<double> X;
- mematrix<int> order;
+ int nids;
+ int ncov;
+ int ngpreds;
+ mematrix<double> weights;
+ mematrix<double> stime;
+ mematrix<int> sstat;
+ mematrix<double> offset;
+ mematrix<int> strata;
+ mematrix<double> X;
+ mematrix<int> order;
- coxph_data(phedata &phed, gendata &gend, int snpnum)
+ coxph_data(phedata &phed, gendata &gend, int snpnum)
+ {
+ nids = gend.nids;
+ ngpreds = gend.ngpreds;
+ if (snpnum>=0)
+ ncov = phed.ncov + ngpreds;
+ else
+ ncov = phed.ncov;
+ if (phed.noutcomes != 2)
{
- nids = gend.nids;
- ngpreds = gend.ngpreds;
- if (snpnum>=0)
- ncov = phed.ncov + ngpreds;
- else
- ncov = phed.ncov;
- if (phed.noutcomes != 2)
- {
- fprintf(stderr,"coxph_data: number of outcomes should be 2 (now: %d)\n",phed.noutcomes);
- exit(1);
- }
-// X.reinit(nids,(ncov+1));
- X.reinit(nids,ncov);
- stime.reinit(nids,1);
- sstat.reinit(nids,1);
- weights.reinit(nids,1);
- offset.reinit(nids,1);
- strata.reinit(nids,1);
- order.reinit(nids,1);
- for (int i=0;i<nids;i++)
- {
+ fprintf(stderr,"coxph_data: number of outcomes should be 2 (now: %d)\n",phed.noutcomes);
+ exit(1);
+ }
+// X.reinit(nids,(ncov+1));
+ X.reinit(nids,ncov);
+ stime.reinit(nids,1);
+ sstat.reinit(nids,1);
+ weights.reinit(nids,1);
+ offset.reinit(nids,1);
+ strata.reinit(nids,1);
+ order.reinit(nids,1);
+ for (int i=0;i<nids;i++)
+ {
// X.put(1.,i,0);
- stime[i] = (phed.Y).get(i,0);
- sstat[i] = int((phed.Y).get(i,1));
- if (sstat[i] != 1 & sstat[i]!=0)
- {
- fprintf(stderr,"coxph_data: status not 0/1 (right order: id, fuptime, status ...)\n",phed.noutcomes);
- exit(1);
- }
- }
- for (int j=0;j<phed.ncov;j++)
- for (int i=0;i<nids;i++)
- X.put((phed.X).get(i,j),i,j);
+ stime[i] = (phed.Y).get(i,0);
+ sstat[i] = int((phed.Y).get(i,1));
+ if (sstat[i] != 1 & sstat[i]!=0)
+ {
+ fprintf(stderr,"coxph_data: status not 0/1 (right order: id, fuptime, status ...)\n",phed.noutcomes);
+ exit(1);
+ }
+ }
+ for (int j=0;j<phed.ncov;j++)
+ for (int i=0;i<nids;i++)
+ X.put((phed.X).get(i,j),i,j);
- if (snpnum>0)
- for (int i=0;i<nids;i++)
- for (int j=0;j<ngpreds;j++)
- X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+j));
+ if (snpnum>0)
+ for (int i=0;i<nids;i++)
+ for (int j=0;j<ngpreds;j++)
+ X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+j));
- for (int i=0;i<nids;i++)
- {
- weights[i] = 1.0;
- offset[i] = 0.0;
- strata[i] = 0;
- }
+ for (int i=0;i<nids;i++)
+ {
+ weights[i] = 1.0;
+ offset[i] = 0.0;
+ strata[i] = 0;
+ }
// sort by time
- double tmptime[nids];
- int passed_sorted[nids];
- for (int i=0;i<nids;i++) {tmptime[i] = stime[i];passed_sorted[i]=0;}
- qsort(tmptime,nids,sizeof(double),cmpfun);
- for (int i=0;i<nids;i++)
- {
- int passed = 0;
- for (int j=0;j<nids;j++)
- if (tmptime[j] == stime[i])
- if (!passed_sorted[j])
- {
- order[i] = j;
- passed_sorted[j] = 1;
- passed = 1;
- break;
- }
- if (passed != 1)
- {
- fprintf(stderr,"can not recover element %d\n",i);
- exit(1);
- }
- }
- stime = reorder(stime,order);
- sstat = reorder(sstat,order);
- weights = reorder(weights,order);
- strata = reorder(strata,order);
- offset = reorder(offset,order);
- X = reorder(X,order);
- X = transpose(X);
+ double tmptime[nids];
+ int passed_sorted[nids];
+ for (int i=0;i<nids;i++) {tmptime[i] = stime[i];passed_sorted[i]=0;}
+ qsort(tmptime,nids,sizeof(double),cmpfun);
+ for (int i=0;i<nids;i++)
+ {
+ int passed = 0;
+ for (int j=0;j<nids;j++)
+ if (tmptime[j] == stime[i])
+ if (!passed_sorted[j])
+ {
+ order[i] = j;
+ passed_sorted[j] = 1;
+ passed = 1;
+ break;
+ }
+ if (passed != 1)
+ {
+ fprintf(stderr,"can not recover element %d\n",i);
+ exit(1);
+ }
+ }
+ stime = reorder(stime,order);
+ sstat = reorder(sstat,order);
+ weights = reorder(weights,order);
+ strata = reorder(strata,order);
+ offset = reorder(offset,order);
+ X = reorder(X,order);
+ X = transpose(X);
// X.print();
// offset.print();
// weights.print();
// stime.print();
// sstat.print();
- }
- void update_snp(gendata &gend, int snpnum)
- {
- // note this sorts by "order"!!!
- for (int i=0;i<nids;i++)
- for (int j=0;j<ngpreds;j++)
- X.put((gend.G).get(i,(snpnum*ngpreds+j)),(ncov-ngpreds+j),order[i]);
- }
- ~coxph_data()
- {
+ }
+ void update_snp(gendata &gend, int snpnum)
+ {
+ // note this sorts by "order"!!!
+ for (int i=0;i<nids;i++)
+ for (int j=0;j<ngpreds;j++)
+ X.put((gend.G).get(i,(snpnum*ngpreds+j)),(ncov-ngpreds+j),order[i]);
+ }
+ ~coxph_data()
+ {
// delete X;
// delete sstat;
// delete stime;
@@ -505,7 +505,7 @@
// delete offset;
// delete strata;
// delete order;
- }
+ }
};
class mlinfo
@@ -524,7 +524,7 @@
mlinfo(char * filename, char * mapname)
{
FILE * infile = fopen(filename,"r");
- if (infile == NULL)
+ if (infile == NULL)
{
fprintf(stderr,"can not open file %s",filename);
exit(1);
@@ -556,7 +556,7 @@
exit(1);
}
for (int i =0;i<7;i++) fscanf(infile,"%s",&tmp);
- for (int i =0;i<nsnps;i++)
+ for (int i =0;i<nsnps;i++)
{
fscanf(infile,"%s",&tmp);
name[i] = tmp;
@@ -580,7 +580,7 @@
std::ifstream instr(mapname);
int BFS = 1000;
char line [BFS], tmp[BFS];
- if (!instr.is_open())
+ if (!instr.is_open())
{
fprintf(stderr,"can not open file %s",filename);
exit(1);
@@ -629,19 +629,19 @@
matrix.reinit(npeople,npeople);
-
+
//idnames[k], if (allmeasured[i]==1)
if (myfile.is_open())
{
- while(myfile.getline(line,MAXIMUM_PEOPLE_AMOUNT))
+ while(myfile.getline(line,MAXIMUM_PEOPLE_AMOUNT))
{
-
-
+
+
std::stringstream line_stream(line);
line_stream >> id;
-
+
if(phe->idnames[row] != id) {std::cerr<<"error:in row "<<row<<" id="<<phe->idnames[row]<<" in inverce variance matrix but id="<<id<<" must be there. Wrong inverce variance matrix (only measured id must be there)\n";exit(1);}
while (line_stream >> val)
@@ -649,7 +649,7 @@
matrix.put(val, row, col);
col++;
}
-
+
if(col != npeople)
{
fprintf(stderr,"error: inv file: Number of columns in row %d equals to %d but people amount is %d\n", row, col, npeople);
@@ -667,7 +667,7 @@
}
- };
+ };
~InvSigma()
@@ -686,14 +686,8 @@
static const unsigned MAXIMUM_PEOPLE_AMOUNT = 1000000;
unsigned npeople; //amount of people
- std::string filename;
+ std::string filename;
mematrix<double> matrix; //file is stored here
};
//________________________________________Maksim_end
-
-
-
-
-
-
Modified: tags/ProbABEL/v.0.1-3/src/reg1.h
===================================================================
--- tags/ProbABEL/v.0.1-3/src/reg1.h 2012-06-11 07:45:12 UTC (rev 915)
+++ tags/ProbABEL/v.0.1-3/src/reg1.h 2012-06-11 20:48:47 UTC (rev 916)
@@ -10,7 +10,7 @@
// last modification: 11-Jan-2009
//
// Author: Yurii S. Aulchenko
-// modified by: Maksim V. Struchalin, 11-Jan-2009
+// modified by: Maksim V. Struchalin, 11-Jan-2009
//
// Modified by Han Chen (hanchen at bu.edu) on Nov 9, 2009
// based on src/reg1.h version 0.2 as of Oct 19, 2009
@@ -29,68 +29,68 @@
// model 3 = recessive 1 df
// model 4 = over-dominant 1 df
{
- if (model==0)
- {
- if(interaction !=0 && !nullmodel)
+ if (model==0)
+ {
+ if(interaction !=0 && !nullmodel)
+ {
+ if(ngpreds == 2)
+ {
+ mematrix<double> nX;
+ nX.reinit(X.nrow,X.ncol+2);
+ int csnp_p1=nX.ncol-4;
+ int csnp_p2=nX.ncol-3;
+ int c1 = nX.ncol-2;
+ int c2 = nX.ncol-1;
+ for (int i=0;i<X.nrow;i++)
+ for (int j=0;j<(X.ncol);j++)
+ nX[i*nX.ncol+j] = X[i*X.ncol+j];
+ for (int i=0;i<nX.nrow;i++)
+ {
+ if (iscox)
+ {
+ nX[i*nX.ncol+c1] = X[i*X.ncol+csnp_p1] * X[i*X.ncol + interaction-1];//Maksim: interaction with SNP;;
+ nX[i*nX.ncol+c2] = X[i*X.ncol+csnp_p2] * X[i*X.ncol + interaction-1];//Maksim: interaction with SNP;;
+ }
+ else
+ {
+ nX[i*nX.ncol+c1] = X[i*X.ncol+csnp_p1] * X[i*X.ncol + interaction];//Maksim: interaction with SNP;;
+ nX[i*nX.ncol+c2] = X[i*X.ncol+csnp_p2] * X[i*X.ncol + interaction];//Maksim: interaction with SNP;;
+ }
+ }
+ //________________________
+ int col_new;
+ if(is_interaction_excluded)
+ {
+ mematrix<double> nX_without_interact_phe;
+ nX_without_interact_phe.reinit(nX.nrow,nX.ncol-1);
+ for(int row=0 ; row < nX.nrow ; row++)
+ {
+//Han Chen
+ col_new=-1;
+ for(int col=0 ; col<nX.ncol ; col++)
{
- if(ngpreds == 2)
- {
- mematrix<double> nX;
- nX.reinit(X.nrow,X.ncol+2);
- int csnp_p1=nX.ncol-4;
- int csnp_p2=nX.ncol-3;
- int c1 = nX.ncol-2;
- int c2 = nX.ncol-1;
- for (int i=0;i<X.nrow;i++)
- for (int j=0;j<(X.ncol);j++)
- nX[i*nX.ncol+j] = X[i*X.ncol+j];
- for (int i=0;i<nX.nrow;i++)
- {
- if (iscox)
- {
- nX[i*nX.ncol+c1] = X[i*X.ncol+csnp_p1] * X[i*X.ncol + interaction-1];//Maksim: interaction with SNP;;
- nX[i*nX.ncol+c2] = X[i*X.ncol+csnp_p2] * X[i*X.ncol + interaction-1];//Maksim: interaction with SNP;;
- }
- else
- {
- nX[i*nX.ncol+c1] = X[i*X.ncol+csnp_p1] * X[i*X.ncol + interaction];//Maksim: interaction with SNP;;
- nX[i*nX.ncol+c2] = X[i*X.ncol+csnp_p2] * X[i*X.ncol + interaction];//Maksim: interaction with SNP;;
- }
- }
- //________________________
- int col_new;
- if(is_interaction_excluded)
- {
- mematrix<double> nX_without_interact_phe;
- nX_without_interact_phe.reinit(nX.nrow,nX.ncol-1);
- for(int row=0 ; row < nX.nrow ; row++)
- {
-//Han Chen
- col_new=-1;
- for(int col=0 ; col<nX.ncol ; col++)
- {
- if(col != interaction && !iscox)
- {
- col_new++;
- nX_without_interact_phe[row*nX_without_interact_phe.ncol+col_new] = nX[row*nX.ncol+col];
- }
- if(col != interaction-1 && iscox)
- {
- col_new++;
- nX_without_interact_phe[row*nX_without_interact_phe.ncol+col_new] = nX[row*nX.ncol+col];
- }
- }//interaction_only, model==0, ngpreds==2
+ if(col != interaction && !iscox)
+ {
+ col_new++;
+ nX_without_interact_phe[row*nX_without_interact_phe.ncol+col_new] = nX[row*nX.ncol+col];
+ }
+ if(col != interaction-1 && iscox)
+ {
+ col_new++;
+ nX_without_interact_phe[row*nX_without_interact_phe.ncol+col_new] = nX[row*nX.ncol+col];
+ }
+ }//interaction_only, model==0, ngpreds==2
//Oct 26, 2009
- }
- return nX_without_interact_phe;
- }//end of is_interaction_excluded
- //________________________
+ }
+ return nX_without_interact_phe;
+ }//end of is_interaction_excluded
+ //________________________
- return(nX);
- }
- if(ngpreds == 1)
+ return(nX);
+ }
+ if(ngpreds == 1)
{
mematrix<double> nX;
nX.reinit(X.nrow,X.ncol+1);
@@ -110,8 +110,8 @@
nX[i*nX.ncol+c1] = X[i*X.ncol+csnp_p1] * X[i*X.ncol + interaction];//Maksim: interaction with SNP;;
}
}
-
+
//________________________
int col_new=-1;
if(is_interaction_excluded)
@@ -123,7 +123,7 @@
col_new=-1;
for(int col=0 ; col<nX.ncol ; col++)
{
- if(col != interaction && !iscox)
+ if(col != interaction && !iscox)
{
col_new++;
nX_without_interact_phe[row*nX_without_interact_phe.ncol+col_new] = nX[row*nX.ncol+col];
@@ -148,7 +148,7 @@
{
return(X);
}
-
+
}
mematrix<double> nX;
if(interaction != 0)
@@ -174,7 +174,7 @@
nX[i*nX.ncol+c1] = X[i*X.ncol+c2];
else if (model==4)
nX[i*nX.ncol+c1] = X[i*X.ncol+c1];
- if(interaction != 0)
+ if(interaction != 0)
nX[i*nX.ncol+c2] = X[i*nX.ncol+interaction] * nX[i*nX.ncol+c1];//Maksim: interaction with SNP
}
//Han Chen
@@ -188,7 +188,7 @@
col_new=-1;
for(int col=0 ; col<nX.ncol ; col++)
{
- if(col != interaction && !iscox)
+ if(col != interaction && !iscox)
{
col_new++;
nX_without_interact_phe[row*nX_without_interact_phe.ncol+col_new] = nX[row*nX.ncol+col];
@@ -261,11 +261,11 @@
//// double sumgt = column_sum(G).get(0,0);
// double meangt = G.column_mean(0);
// double meany = Y.column_mean(0);
-//
+//
// G = G - meangt;
// Y = Y - meany;
// double U, V;
-// mematrix<double> U_mema, V_mema;
+// mematrix<double> U_mema, V_mema;
//
//
// U_mema = transpose(G)*invvarmatrix;
@@ -284,11 +284,11 @@
// double sebe=sqrt(1./V);
// sebeta.put(sebe,1,0);
//
-//
//
//
+//
// /*
-// for(snp in snps)
+// for(snp in snps)
// {
// print(snp)
// ns <- which(mm$snpnames==snp)
@@ -301,12 +301,12 @@
// chi2 <- UU*UU/VV
// print(c(re$chi2.1df[ns],qt$chi2.1df[ns],gr$chi2.1df[ns],mm$chi2.1df[ns],(beta/se)^2))
// print(c(re$effB[ns],qt$effB[ns],gr$effB[ns],mm$effB[ns],beta))
-// }
+// }
//
// */
-//
//
//
+//
///*
// int length_beta = (rdata.X).ncol;
// beta.reinit(length_beta,1);
@@ -339,35 +339,35 @@
//suda ineraction parameter
// model should come here
mematrix<double> X = apply_model(rdata.X,model, interaction, ngpreds, false, nullmodel);
-
+
int length_beta = X.ncol;
beta.reinit(length_beta,1);
sebeta.reinit(length_beta,1);
//Han Chen
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);}
- }
+ {covariance.reinit(length_beta-2,1);}
+ else
+ {covariance.reinit(length_beta-1,1);}
+ }
//Oct 26, 2009
mematrix<double> tX = transpose(X);
-
- if(invvarmatrix.nrow!=0 && invvarmatrix.ncol!=0)
+
+ if(invvarmatrix.nrow!=0 && invvarmatrix.ncol!=0)
{
- tX = tX*invvarmatrix;
+ tX = tX*invvarmatrix;
// X = invvarmatrix*X; std::cout<<"new tX.nrow="<<X.nrow<<" tX.ncol="<<X.ncol<<"\n";
}
-
-
-
+
+
+
mematrix<double> tXX = tX*X;
// double N = tXX.get(0,0);
double N = X.nrow;
if (verbose) {printf("tXX:\n");tXX.print();}
- //
+ //
// use cholesky to invert
//
mematrix<double> tXX_i = tXX;
@@ -408,19 +408,19 @@
// now compute residual variance
// sigma2 = 0.;
-// for (int i =0;i<(rdata.Y).nrow;i++)
+// for (int i =0;i<(rdata.Y).nrow;i++)
// sigma2 += ((rdata.Y).get(i,0))*((rdata.Y).get(i,0));
-// for (int i=0;i<length_beta;i++)
+// for (int i=0;i<length_beta;i++)
// sigma2 -= 2. * (beta.get(i,0)) * tXY.get(i,0);
-// for (int i=0;i<(length_beta);i++)
-// for (int j=0;j<(length_beta);j++)
+// for (int i=0;i<(length_beta);i++)
+// for (int j=0;j<(length_beta);j++)
// sigma2 += (beta.get(i,0)) * (beta.get(j,0)) * tXX.get(i,j);
// std::cout<<"sigma2="<<sigma2<<"\n";
// std::cout<<"sigma2_internal="<<sigma2_internal<<"\n";
// replaced for ML
-// sigma2_internal = sigma2/(N - double(length_beta) - 1);
+// sigma2_internal = sigma2/(N - double(length_beta) - 1);
sigma2 /= N;
// std::cout<<"N="<<N<<", length_beta="<<length_beta<<"\n";
@@ -468,24 +468,24 @@
}
for (int i=0;i<(length_beta);i++)
- {
+ {
if (robust) {
double value = sqrt(robust_sigma2.get(i,i));
sebeta.put(value,i,0);
//Han Chen
if (i>0)
{
- if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
- {if (i>1)
- {double covval=robust_sigma2.get(i,i-2);
- covariance.put(covval,i-2,0);}
- }
- else
- {
- double covval=robust_sigma2.get(i,i-1);
- covariance.put(covval,i-1,0);
- }
- }
+ if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
+ {if (i>1)
+ {double covval=robust_sigma2.get(i,i-2);
+ covariance.put(covval,i-2,0);}
+ }
+ else
+ {
+ double covval=robust_sigma2.get(i,i-1);
+ covariance.put(covval,i-1,0);
+ }
+ }
//Oct 26, 2009
} else {
double value = sqrt(sigma2_internal*tXX_i.get(i,i));
@@ -493,17 +493,17 @@
//Han Chen
if (i>0)
{
- if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
- {if (i>1)
- {double covval=sigma2_internal*tXX_i.get(i,i-2);
- covariance.put(covval,i-2,0);}
- }
- else
- {
- double covval=sigma2_internal*tXX_i.get(i,i-1);
- covariance.put(covval,i-1,0);
- }
- }
+ if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
+ {if (i>1)
+ {double covval=sigma2_internal*tXX_i.get(i,i-2);
+ covariance.put(covval,i-2,0);}
+ }
+ else
+ {
+ double covval=sigma2_internal*tXX_i.get(i,i-1);
+ covariance.put(covval,i-1,0);
+ }
+ }
//Oct 26, 2009
}
}
@@ -516,13 +516,13 @@
beta.reinit(X.ncol,1);
sebeta.reinit(X.ncol,1);
double N = double(resid.nrow);
-
+
mematrix<double> tX = transpose(X);
-
+
if(invvarmatrix.nrow!=0 && invvarmatrix.ncol!=0) tX = tX*invvarmatrix;
//if(invvarmatrix.nrow!=0 && invvarmatrix.ncol!=0) X = invvarmatrix*X;
-
-
+
+
mematrix<double> u = tX*resid;
mematrix<double> v = tX*X;
mematrix<double> csum = column_sum(X);
@@ -538,7 +538,7 @@
beta = v_i*u;
double sr = 0.;
double srr =0.;
- for (int i=0;i<resid.nrow;i++)
+ for (int i=0;i<resid.nrow;i++)
{
sr += resid[i];
srr += resid[i]*resid[i];
@@ -598,15 +598,15 @@
//Han Chen
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);}
- }
+ {covariance.reinit(length_beta-2,1);}
+ else
+ {covariance.reinit(length_beta-1,1);}
+ }
//Oct 26, 2009
mematrix<double> W((X).nrow,1);
mematrix<double> z((X).nrow,1);
mematrix<double> tXWX(length_beta,length_beta);
- mematrix<double> tXWX_i(length_beta,length_beta);
+ mematrix<double> tXWX_i(length_beta,length_beta);
mematrix<double> tXWz(length_beta,1);
double prev = (rdata.Y).column_mean(0);
@@ -652,7 +652,7 @@
N = tXWX.get(0,0);
if (verbose) {printf("tXWX:\n");tXWX.print();}
- //
+ //
// use cholesky to invert
//
tXWX_i = tXWX;
@@ -690,44 +690,44 @@
}
for (int i=0;i<(length_beta);i++)
- {
- if (robust)
- {
+ {
+ if (robust)
+ {
double value = sqrt(robust_sigma2.get(i,i));
sebeta.put(value,i,0);
//Han Chen
if (i>0)
{
- if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
- {if (i>1)
- {double covval=robust_sigma2.get(i,i-2);
- covariance.put(covval,i-2,0);}
- }
- else
- {
- double covval=robust_sigma2.get(i,i-1);
- covariance.put(covval,i-1,0);
- }
- }
+ if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
+ {if (i>1)
+ {double covval=robust_sigma2.get(i,i-2);
+ covariance.put(covval,i-2,0);}
+ }
+ else
+ {
+ double covval=robust_sigma2.get(i,i-1);
+ covariance.put(covval,i-1,0);
+ }
+ }
//Oct 26, 2009
- }
- else {
+ }
+ else {
double value = sqrt(tXWX_i.get(i,i));
sebeta.put(value,i,0);
//Han Chen
if (i>0)
{
- if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
- {if (i>1)
- {double covval=tXWX_i.get(i,i-2);
- covariance.put(covval,i-2,0);}
- }
- else
- {
- double covval=tXWX_i.get(i,i-1);
- covariance.put(covval,i-1,0);
- }
- }
+ if (model==0 && interaction!=0 && ngpreds==2 && length_beta>2)
+ {if (i>1)
+ {double covval=tXWX_i.get(i,i-2);
+ covariance.put(covval,i-2,0);}
+ }
+ else
+ {
+ double covval=tXWX_i.get(i,i-1);
+ covariance.put(covval,i-1,0);
+ }
+ }
//Oct 26, 2009
}
}
@@ -741,10 +741,10 @@
beta.reinit(X.ncol,1);
sebeta.reinit(X.ncol,1);
double N = double(resid.nrow);
-
- mematrix<double> tX=transpose(X);
+
+ mematrix<double> tX=transpose(X);
if(invvarmatrix.nrow!=0 && invvarmatrix.ncol!=0) tX = tX*invvarmatrix;
-
+
mematrix<double> u = tX*resid;
mematrix<double> v = tX*X;
mematrix<double> csum = column_sum(X);
@@ -760,7 +760,7 @@
beta = v_i*u;
double sr = 0.;
double srr =0.;
- for (int i=0;i<resid.nrow;i++)
+ for (int i=0;i<resid.nrow;i++)
{
sr += resid[i];
srr += resid[i]*resid[i];
@@ -776,68 +776,68 @@
};
-void coxfit2(int *maxiter, int *nusedx, int *nvarx,
- double *time, int *status, double *covar2,
+void coxfit2(int *maxiter, int *nusedx, int *nvarx,
+ double *time, int *status, double *covar2,
double *offset, double *weights, int *strata,
- double *means, double *beta, double *u,
- double *imat2, double loglik[2], int *flag,
+ double *means, double *beta, double *u,
+ double *imat2, double loglik[2], int *flag,
double *work, double *eps, double *tol_chol,
double *sctest);
class coxph_reg
{
public:
- mematrix<double> beta;
- mematrix<double> sebeta;
- mematrix<double> residuals;
- double sigma2;
- double loglik;
- double chi2_score;
- int niter;
+ mematrix<double> beta;
+ mematrix<double> sebeta;
+ mematrix<double> residuals;
+ double sigma2;
+ double loglik;
+ double chi2_score;
+ int niter;
- coxph_reg(coxph_data &cdata)
- {
- beta.reinit(cdata.X.nrow,1);
- sebeta.reinit(cdata.X.nrow,1);
- loglik=-9.999e+32;
- sigma2=-1.;
- chi2_score=-1.;
- niter=0;
- }
- ~coxph_reg()
- {
+ coxph_reg(coxph_data &cdata)
+ {
+ beta.reinit(cdata.X.nrow,1);
+ sebeta.reinit(cdata.X.nrow,1);
+ loglik=-9.999e+32;
+ sigma2=-1.;
+ chi2_score=-1.;
+ niter=0;
+ }
+ ~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)
+ {
// cout << "model = " << model << "\n";
// cdata.X.print();
- mematrix<double> X = t_apply_model(cdata.X,model, interaction, ngpreds, iscox, nullmodel);
+ mematrix<double> X = t_apply_model(cdata.X,model, interaction, ngpreds, iscox, nullmodel);
// X.print();
- int length_beta = X.nrow;
- beta.reinit(length_beta,1);
- sebeta.reinit(length_beta,1);
- mematrix<double> newoffset = cdata.offset;
- newoffset = cdata.offset - (cdata.offset).column_mean(0);
- mematrix<double> means(X.nrow,1);
- for (int i=0;i<X.nrow;i++) beta[i]=0.;
- mematrix<double> u(X.nrow,1);
- mematrix<double> imat(X.nrow,X.nrow);
- double work[X.ncol*2+2*(X.nrow)*(X.nrow)+3*(X.nrow)];
- double loglik_int[2];
- int flag;
- double sctest=1.0;
-
- coxfit2(&maxiter,&cdata.nids,&X.nrow,
- cdata.stime.data,cdata.sstat.data,X.data,
- newoffset.data,cdata.weights.data,cdata.strata.data,
- means.data,beta.data,u.data,
- 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));
- loglik = loglik_int[1];
- niter = maxiter;
- }
+ int length_beta = X.nrow;
+ beta.reinit(length_beta,1);
+ sebeta.reinit(length_beta,1);
+ mematrix<double> newoffset = cdata.offset;
+ newoffset = cdata.offset - (cdata.offset).column_mean(0);
+ mematrix<double> means(X.nrow,1);
+ for (int i=0;i<X.nrow;i++) beta[i]=0.;
+ mematrix<double> u(X.nrow,1);
+ mematrix<double> imat(X.nrow,X.nrow);
+ double work[X.ncol*2+2*(X.nrow)*(X.nrow)+3*(X.nrow)];
+ double loglik_int[2];
+ int flag;
+ double sctest=1.0;
+
+ coxfit2(&maxiter,&cdata.nids,&X.nrow,
+ cdata.stime.data,cdata.sstat.data,X.data,
+ newoffset.data,cdata.weights.data,cdata.strata.data,
+ means.data,beta.data,u.data,
+ 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));
+ loglik = loglik_int[1];
+ niter = maxiter;
+ }
};
More information about the Genabel-commits
mailing list