Summary: 0.99 in progress

diff --git a/src/imach.c b/src/imach.c
index 76ce1ef10456338b6d581ada89b1fc2668709e49..f67de2a4fddb883b65c30010030514702ac28715 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.252  2016/09/15 21:15:37  brouard
+  *** empty log message ***
+
   Revision 1.251  2016/09/15 15:01:13  brouard
   Summary: not working
 
@@ -4186,17 +4189,84 @@ void pstamp(FILE *fichier)
   fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
 }
 
+int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) {
+
+  /* y=a+bx regression */
+  double   sumx = 0.0;                        /* sum of x                      */
+  double   sumx2 = 0.0;                       /* sum of x**2                   */
+  double   sumxy = 0.0;                       /* sum of x * y                  */
+  double   sumy = 0.0;                        /* sum of y                      */
+  double   sumy2 = 0.0;                       /* sum of y**2                   */
+  double   sume2; /* sum of square or residuals */
+  double yhat;
+  
+  double denom=0;
+  int i;
+  int ne=*no;
+  
+  for ( i=ifi, ne=0;i<=ila;i++) {
+    if(!isfinite(x[i]) || !isfinite(y[i])){
+      /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+      continue;
+    }
+    ne=ne+1;
+    sumx  += x[i];       
+    sumx2 += x[i]*x[i];  
+    sumxy += x[i] * y[i];
+    sumy  += y[i];      
+    sumy2 += y[i]*y[i]; 
+    denom = (ne * sumx2 - sumx*sumx);
+    /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+  } 
+  
+  denom = (ne * sumx2 - sumx*sumx);
+  if (denom == 0) {
+    // vertical, slope m is infinity
+    *b = INFINITY;
+    *a = 0;
+    if (r) *r = 0;
+    return 1;
+  }
+  
+  *b = (ne * sumxy  -  sumx * sumy) / denom;
+  *a = (sumy * sumx2  -  sumx * sumxy) / denom;
+  if (r!=NULL) {
+    *r = (sumxy - sumx * sumy / ne) /          /* compute correlation coeff     */
+      sqrt((sumx2 - sumx*sumx/ne) *
+          (sumy2 - sumy*sumy/ne));
+  }
+  *no=ne;
+  for ( i=ifi, ne=0;i<=ila;i++) {
+    if(!isfinite(x[i]) || !isfinite(y[i])){
+      /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+      continue;
+    }
+    ne=ne+1;
+    yhat = y[i] - *a -*b* x[i];
+    sume2  += yhat * yhat ;       
+    
+    denom = (ne * sumx2 - sumx*sumx);
+    /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+  } 
+  *sb = sqrt(sume2/(ne-2)/(sumx2 - sumx * sumx /ne));
+  *sa= *sb * sqrt(sumx2/ne);
+  
+  return 0; 
+}
+
 /************ Frequencies ********************/
 void  freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
                  int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \
                  int firstpass,  int lastpass, int stepm, int weightopt, char model[])
 {  /* Some frequencies as well as proposing some starting values */
   
-  int i, m, jk, j1, bool, z1,j, k, iv, jj=0;
+  int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0;
   int iind=0, iage=0;
   int mi; /* Effective wave */
   int first;
   double ***freq; /* Frequencies */
+  double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */
+  int no;
   double *meanq;
   double **meanqt;
   double *pp, **prop, *posprop, *pospropt;
@@ -4249,6 +4319,8 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
   }
   fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
   
+  y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
+  x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+4+AGEMARGE);
   j1=0;
   
@@ -4266,7 +4338,15 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
   dateintsum=0;
   k2cpt=0;
 
-  for (j = 0; j <= cptcoveff; j+=cptcoveff){   /* j= 0 constant model */
+  if(cptcoveff == 0 )
+    nl=1;  /* Constant model only */
+  else
+    nl=2;
+  for (nj = 1; nj <= nl; nj++){   /* nj= 1 constant model, nl number of loops. */
+    if(nj==1)
+      j=0;  /* First pass for the constant */
+    else
+      j=cptcoveff; /* Other passes for the covariate values */
     first=1;
     for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */
       posproptt=0.;
@@ -4512,6 +4592,7 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
              if(first==1){
                printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
              }
+             /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); */
              fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
            }
            if(jk!=0 && m!=0)
@@ -4559,7 +4640,7 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
          for(k=1; k <=(nlstate+ndeath); k++){
            if (k != i) {
              for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */
-               if(jj==1){  /* Constant case */
+               if(jj==1){  /* Constant case (in fact cste + age) */
                  if(j1==1){ /* All dummy covariates to zero */
                    freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */
                    freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */
@@ -4569,8 +4650,15 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
                    fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
                    pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
                  }
-               }else if(jj==2 || nagesqr==1){ /* age or age*age parameter */
-                 ;
+               }else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */
+                 for(iage=iagemin; iage <= iagemax+3; iage++){
+                   x[iage]= (double)iage;
+                   y[iage]= log(freq[i][k][iage]/freq[i][i][iage]);
+                   /* printf("i=%d, k=%d, jk=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,jk,j1,jj, iage, y[iage]); */
+                 }
+                 linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */
+                 pstart[jk]=b;
+                 pstart[jk-1]=a;
                }else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj)  && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ 
                  printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
                  printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
@@ -4599,7 +4687,9 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
            printf("%d%d ",i,k);
            fprintf(ficlog,"%d%d ",i,k);
            for(jj=1; jj <=ncovmodel; jj++){
-             if(jj==1){
+             pstart[jk]=p[jk]; /* Setting pstart to p values by default */
+             if(jj==1){ /* Age has to be done */
+               pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
                printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
                fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
              }
@@ -4631,7 +4721,7 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
     } /* end j=0 */
   } /* end j */
 
-  if(mle == -2){
+  if(mle == -2){  /* We want to use these values as starting values */
     for(i=1, jk=1; i <=nlstate; i++){
       for(j=1; j <=nlstate+ndeath; j++){
        if(j!=i){
@@ -4659,6 +4749,8 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
   fclose(ficresphtmfr);
   free_vector(meanq,1,nqfveff);
   free_matrix(meanqt,1,lastpass,1,nqtveff);
+  free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
+  free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_vector(pospropt,1,nlstate);
   free_vector(posprop,1,nlstate);
@@ -6310,7 +6402,7 @@ void printinghtml(char fileresu[], char title[], char datafile[], int firstpass,
 
    for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
-     if(TKresult[nres]!= k1)
+     if(m != 1 && TKresult[nres]!= k1)
        continue;
 
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
@@ -6442,7 +6534,7 @@ See page 'Matrix of variance-covariance of one-step probabilities' below. \n", r
 
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
    for(k1=1; k1<=m;k1++){
-     if(TKresult[nres]!= k1)
+     if(m != 1 && TKresult[nres]!= k1)
        continue;
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
@@ -6538,7 +6630,7 @@ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar,
     for (k1=1; k1<= m ; k1 ++){ /* For each valid combination of covariate */
       for(nres=1; nres <= nresult; nres++){ /* For each resultline */
        /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-       if(TKresult[nres]!= k1)
+       if(m != 1 && TKresult[nres]!= k1)
          continue;
        /* We are interested in selected combination by the resultline */
        /* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */
@@ -6620,7 +6712,7 @@ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar,
   /*2 eme*/
   for (k1=1; k1<= m ; k1 ++){  
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
        continue;
       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
@@ -6680,7 +6772,7 @@ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar,
   /*3eme*/
   for (k1=1; k1<= m ; k1 ++){
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
        continue;
 
       for (cpt=1; cpt<= nlstate ; cpt ++) {
@@ -6729,7 +6821,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subd
   /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
        continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
        fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
@@ -6775,7 +6867,7 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1++){ /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
        continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
        fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
@@ -6829,7 +6921,7 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* CV preval stable (period) for each covariate */
   for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    if(TKresult[nres]!= k1)
+    if(m != 1 && TKresult[nres]!= k1)
       continue;
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       
@@ -6876,7 +6968,7 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
     /* CV back preval stable (period) for each covariate */
     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
        continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
        fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
@@ -6927,7 +7019,7 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
     
     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
        continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
        fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
@@ -7058,7 +7150,7 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
     fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
     for(jk=1; jk <=m; jk++)  /* For each combination of covariate */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= jk)
+      if(m != 1 && TKresult[nres]!= jk)
        continue;
       fprintf(ficgp,"# Combination of dummy  jk=%d and ",jk);
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
@@ -7420,7 +7512,7 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
 /*           if (h==(int)(YEARM*yearp)){ */
   for(nres=1; nres <= nresult; nres++) /* For each resultline */
   for(k=1; k<=i1;k++){
-    if(TKresult[nres]!= k)
+    if(i1 != 1 && TKresult[nres]!= k)
       continue;
     if(invalidvarcomb[k]){
       printf("\nCombination (%d) projection ignored because no cases \n",k); 
@@ -8371,6 +8463,11 @@ int decoderesult ( char resultline[], int nres)
   if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */
   }
+  if(j == 0){ /* Resultline but no = */
+    TKresult[nres]=0; /* Combination for the nresult and the model */
+    return (0);
+  }
+    
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
     printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
     fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
@@ -9307,7 +9404,7 @@ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxp
 
   for(k=1; k<=i1;k++){ /* For each combination k of dummy covariates in the model */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
        continue;
 
       /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
@@ -9404,7 +9501,7 @@ int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double a
   
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+     if(i1 != 1 && TKresult[nres]!= k)
        continue;
       //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
       fprintf(ficresplb,"#******");
@@ -9507,7 +9604,7 @@ int hPijx(double *p, int bage, int fage){
                /*      k=k+1;  */
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
        continue;
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
@@ -9586,7 +9683,7 @@ int hPijx(double *p, int bage, int fage){
   /*   k=k+1;  */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
        continue;
       fprintf(ficrespijb,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)
@@ -11094,9 +11191,10 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
     while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
       if (num_filled == 0){
        resultline[0]='\0';
+       printf("Warning %d: no result line should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line);
       break;
       } else if (num_filled != 1){
-       printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
+       printf("ERROR %d: result line should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line);
       }
       nresult++; /* Sum of resultlines */
       printf("Result %d: result=%s\n",nresult, resultline);
@@ -11268,7 +11366,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
     
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
        continue;
       fprintf(ficreseij,"\n#****** ");
       printf("\n#****** ");
@@ -11341,7 +11439,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
     
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
        continue;
       printf("\n#****** Result for:");
       fprintf(ficrest,"\n#****** Result for:");
@@ -11482,7 +11580,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
 
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
        continue;
       fprintf(ficresvpl,"\n#****** ");
       printf("\n#****** ");