--- imach/src/imach.c 2015/11/18 17:41:20 1.210 +++ imach/src/imach.c 2015/12/16 08:52:24 1.215 @@ -1,6 +1,23 @@ -/* $Id: imach.c,v 1.210 2015/11/18 17:41:20 brouard Exp $ +/* $Id: imach.c,v 1.215 2015/12/16 08:52:24 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.215 2015/12/16 08:52:24 brouard + Summary: 0.98r4 working + + Revision 1.214 2015/12/16 06:57:54 brouard + Summary: temporary not working + + Revision 1.213 2015/12/11 18:22:17 brouard + Summary: 0.98r4 + + Revision 1.212 2015/11/21 12:47:24 brouard + Summary: minor typo + + Revision 1.211 2015/11/21 12:41:11 brouard + Summary: 0.98r3 with some graph of projected cross-sectional + + Author: Nicolas Brouard + Revision 1.210 2015/11/18 17:41:20 brouard Summary: Start working on projected prevalences @@ -757,6 +774,8 @@ typedef struct { #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 +/*#define decodtabm(h,k,cptcoveff)= (h <= (1<> (k-1)) & 1) +1 : -1)*/ +#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 #define MAXN 20000 #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 @@ -773,12 +792,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.210 2015/11/18 17:41:20 brouard Exp $ */ +/* $Id: imach.c,v 1.215 2015/12/16 08:52:24 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.210 $ $Date: 2015/11/18 17:41:20 $"; +char fullversion[]="$Revision: 1.215 $ $Date: 2015/12/16 08:52:24 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -814,7 +833,7 @@ double **matprod2(); /* test */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ /*FILE *fic ; */ /* Used in readdata only */ -FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; +FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; int globpr=0; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ @@ -1376,7 +1395,30 @@ char *subdirf3(char fileres[], char *pre strcat(tmpout,fileres); return tmpout; } + +/*************** function subdirfext ***********/ +char *subdirfext(char fileres[], char *preop, char *postop) +{ + + strcpy(tmpout,preop); + strcat(tmpout,fileres); + strcat(tmpout,postop); + return tmpout; +} +/*************** function subdirfext3 ***********/ +char *subdirfext3(char fileres[], char *preop, char *postop) +{ + + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); + strcat(tmpout,preop); + strcat(tmpout,fileres); + strcat(tmpout,postop); + return tmpout; +} + char *asc_diff_time(long time_sec, char ascdiff[]) { long sec_left, days, hours, minutes; @@ -2213,6 +2255,7 @@ double ***hpxij(double ***po, int nhstep double **out, cov[NCOVMAX+1]; double **newm; double agexact; + double agebegin, ageend; /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) @@ -2226,7 +2269,7 @@ double ***hpxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; + agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ cov[2]=agexact; if(nagesqr==1) cov[3]= agexact*agexact; @@ -2597,6 +2640,7 @@ double funcone( double *x) int s1, s2; double bbh, survp; double agexact; + double agebegin, ageend; /*extern weight */ /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ @@ -2615,7 +2659,12 @@ double funcone( double *x) oldm[ii][j]=(ii==j ? 1.0 : 0.0); savm[ii][j]=(ii==j ? 1.0 : 0.0); } - for(d=0; d- Probability p%dj by origin %d and destination j %s-p%dj.png
\ + fprintf(fichtm,"
- Probability p%dj by origin %d and destination j. Dot's sizes are related to corresponding weight: %s-p%dj.png
\ ",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); } fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ @@ -3178,25 +3227,61 @@ void pstamp(FILE *fichier) } /************ Frequencies ********************/ -void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) +void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ + int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],\ + int firstpass, int lastpass, int stepm, int weightopt, char model[]) { /* Some frequencies */ int i, m, jk, j1, bool, z1,j; + int mi; /* Effective wave */ int first; double ***freq; /* Frequencies */ double *pp, **prop; double pos,posprop, k2, dateintsum=0,k2cpt=0; - char fileresp[FILENAMELENGTH]; - + char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; + double agebegin, ageend; + pp=vector(1,nlstate); prop=matrix(1,nlstate,iagemin,iagemax+3); strcpy(fileresp,"P_"); strcat(fileresp,fileresu); + /*strcat(fileresphtm,fileresu);*/ if((ficresp=fopen(fileresp,"w"))==NULL) { printf("Problem with prevalence resultfile: %s\n", fileresp); fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } + + strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); + if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { + printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtm,"\nIMaCh PHTM_ %s\n %s
%s \ +
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ + fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies and prevalence by age at begin of transition

\n",fileresphtm, fileresphtm); + + strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); + if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { + printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s \ +
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ + fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions by age at begin of transition

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); j1=0; @@ -3205,10 +3290,7 @@ void freqsummary(char fileres[], int ia first=1; - /* for(k1=1; k1<=j ; k1++){ */ /* Loop on covariates */ - /* for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */ - /* j1++; */ - for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ + for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */ /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ for (i=-5; i<=nlstate+ndeath; i++) @@ -3222,7 +3304,7 @@ void freqsummary(char fileres[], int ia dateintsum=0; k2cpt=0; - for (i=1; i<=imx; i++) { + for (i=1; i<=imx; i++) { /* For each individual i */ bool=1; if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ for (z1=1; z1<=cptcoveff; z1++) @@ -3235,25 +3317,38 @@ void freqsummary(char fileres[], int ia /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ } } /* cptcovn > 0 */ - + if (bool==1){ - for(m=firstpass; m<=lastpass; m++){ - k2=anint[m][i]+(mint[m][i]/12.); - /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ - if(agev[m][i]==0) agev[m][i]=iagemax+1; - if(agev[m][i]==1) agev[m][i]=iagemax+2; - if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; + /* for(m=firstpass; m<=lastpass; m++){ */ + for(mi=1; mi=firstpass && m <=lastpass){ + k2=anint[m][i]+(mint[m][i]/12.); + /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ + if(agev[m][i]==0) agev[m][i]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ + if(agev[m][i]==1) agev[m][i]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ + if (s[m][i]>0 && s[m][i]<=nlstate) /* If status at wave m is known and a live state */ + prop[s[m][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */ if (m1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) { - dateintsum=dateintsum+k2; - k2cpt++; - /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ + /* if(s[m][i]==4 && s[m+1][i]==4) */ + /* printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i]); */ + if(s[m][i]==-1) + printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i],agebegin, ageend, (int)((agebegin+ageend)/2.)); + freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */ + /* freq[s[m][i]][s[m+1][i]][(int)((agebegin+ageend)/2.)] += weight[i]; */ + freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */ } - /*}*/ + } + if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) { + dateintsum=dateintsum+k2; + k2cpt++; + /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ + } + /*}*/ } /* end m */ } /* end bool */ } /* end i = 1 to imx */ @@ -3262,24 +3357,57 @@ void freqsummary(char fileres[], int ia pstamp(ficresp); if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "\n

********** Variable "); + fprintf(ficresphtmfr, "\n

********** Variable "); + for (z1=1; z1<=cptcoveff; z1++){ + fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + } + fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "**********

\n"); + fprintf(ficresphtmfr, "**********\n"); fprintf(ficlog, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficlog, "**********\n#"); + fprintf(ficlog, "**********\n"); } - for(i=1; i<=nlstate;i++) + fprintf(ficresphtm,""); + for(i=1; i<=nlstate;i++) { fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); + fprintf(ficresphtm, "",i,i); + } fprintf(ficresp, "\n"); + fprintf(ficresphtm, "\n"); + /* Header of frequency table by age */ + fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); + fprintf(ficresphtmfr," "); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",jk,m); + } + } + fprintf(ficresphtmfr, "\n"); + + /* For each age */ for(i=iagemin; i <= iagemax+3; i++){ - if(i==iagemax+3){ + fprintf(ficresphtm,""); + if(i==iagemax+1){ + fprintf(ficlog,"1"); + fprintf(ficresphtmfr," "); + }else if(i==iagemax+2){ + fprintf(ficlog,"0"); + fprintf(ficresphtmfr," "); + }else if(i==iagemax+3){ fprintf(ficlog,"Total"); + fprintf(ficresphtmfr," "); }else{ if(first==1){ first=0; printf("See log file for details...\n"); } + fprintf(ficresphtmfr," ",i); fprintf(ficlog,"Age %d", i); } for(jk=1; jk <=nlstate ; jk++){ @@ -3322,32 +3450,47 @@ void freqsummary(char fileres[], int ia if( i <= iagemax){ if(pos>=1.e-5){ fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); + fprintf(ficresphtm,"",i,prop[jk][i]/posprop, prop[jk][i],posprop); /*probs[i][jk][j1]= pp[jk]/pos;*/ /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ } - else + else{ fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); + fprintf(ficresphtm,"",i, prop[jk][i],posprop); + } } } - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) - if(freq[jk][m][i] !=0 ) { - if(first==1) - printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[jk][m][i] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); + } fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]); } - if(i <= iagemax) + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",freq[jk][m][i]); + } + } + fprintf(ficresphtmfr,"\n "); + if(i <= iagemax){ fprintf(ficresp,"\n"); + fprintf(ficresphtm,"\n"); + } if(first==1) printf("Others in log...\n"); fprintf(ficlog,"\n"); } /* end loop i */ + fprintf(ficresphtm,"
Age%d%d
0
Unknown
Total
%d%d%.5f%.0f%.0f%dNaNq%.0f%.0f%.0f
\n"); + fprintf(ficresphtmfr,"\n"); /*}*/ } /* end j1 */ dateintmean=dateintsum/k2cpt; fclose(ficresp); + fclose(ficresphtm); + fclose(ficresphtmfr); free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); free_vector(pp,1,nlstate); free_matrix(prop,1,nlstate,iagemin, iagemax+3); @@ -3363,6 +3506,9 @@ void prevalence(double ***probs, double */ int i, m, jk, j1, bool, z1,j; + int mi; /* Effective wave */ + int iage; + double agebegin, ageend; double **prop; double posprop; @@ -3382,54 +3528,57 @@ void prevalence(double ***probs, double first=1; for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ - /*for(i1=1; i1<=ncodemax[k1];i1++){ - j1++;*/ - - for (i=1; i<=nlstate; i++) - for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0.0; - - for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) - bool=0; - } - if (bool==1) { - for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/ + for (i=1; i<=nlstate; i++) + for(iage=iagemin; iage <= iagemax+3; iage++) + prop[i][iage]=0.0; + + for (i=1; i<=imx; i++) { /* Each individual */ + bool=1; + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcoveff; z1++) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) + bool=0; + } + if (bool==1) { + /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ + for(mi=1; mi=firstpass && m <=lastpass){ y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ if(agev[m][i]==0) agev[m][i]=iagemax+1; if(agev[m][i]==1) agev[m][i]=iagemax+2; if((int)agev[m][i] iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); - if (s[m][i]>0 && s[m][i]<=nlstate) { + if (s[m][i]>0 && s[m][i]<=nlstate) { /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ - prop[s[m][i]][(int)agev[m][i]] += weight[i]; - prop[s[m][i]][iagemax+3] += weight[i]; - } - } + prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ + prop[s[m][i]][iagemax+3] += weight[i]; + } /* end valid statuses */ + } /* end selection of dates */ } /* end selection of waves */ - } - } - for(i=iagemin; i <= iagemax+3; i++){ - for(jk=1,posprop=0; jk <=nlstate ; jk++) { - posprop += prop[jk][i]; - } - - for(jk=1; jk <=nlstate ; jk++){ - if( i <= iagemax){ - if(posprop>=1.e-5){ - probs[i][jk][j1]= prop[jk][i]/posprop; - } else{ - if(first==1){ - first=0; - printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); - } + } /* end effective waves */ + } /* end bool */ + } + for(i=iagemin; i <= iagemax+3; i++){ + for(jk=1,posprop=0; jk <=nlstate ; jk++) { + posprop += prop[jk][i]; + } + + for(jk=1; jk <=nlstate ; jk++){ + if( i <= iagemax){ + if(posprop>=1.e-5){ + probs[i][jk][j1]= prop[jk][i]/posprop; + } else{ + if(first==1){ + first=0; + printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); } - } - }/* end jk */ - }/* end i */ + } + } + }/* end jk */ + }/* end i */ /*} *//* end i1 */ } /* end j1 */ @@ -3452,17 +3601,18 @@ void concatwav(int wav[], int **dh, int int i, mi, m; /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; double sum=0., jmean=0.;*/ - int first; + int first, firstwo; int j, k=0,jk, ju, jl; double sum=0.; first=0; + firstwo=0; jmin=100000; jmax=-1; jmean=0.; - for(i=1; i<=imx; i++){ + for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ mi=0; m=firstpass; - while(s[m][i] <= nlstate){ + while(s[m][i] <= nlstate){ /* a live state */ if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5) mw[++mi][i]=m; if(m >=lastpass) @@ -3470,13 +3620,26 @@ void concatwav(int wav[], int **dh, int else m++; }/* end while */ - if (s[m][i] > nlstate){ + if (s[m][i] > nlstate){ /* In a death state */ mi++; /* Death is another wave */ /* if(mi==0) never been interviewed correctly before death */ /* Only death is a correct wave */ mw[mi][i]=m; + }else if (andc[i] != 9999) { /* A death occured after lastpass */ + m++; + mi++; + s[m][i]=nlstate+1; /* We are setting the status to the last of non live state */ + mw[mi][i]=m; + nbwarn++; + if(firstwo==0){ + printf("Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\nOthers in log file only\n",num[i],i,lastpass,nlstate+1, m); + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m); + firstwo=1; + } + if(firstwo==1){ + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m); + } } - wav[i]=mi; if(mi==0){ nbwarn++; @@ -3489,7 +3652,9 @@ void concatwav(int wav[], int **dh, int } } /* end mi==0 */ } /* End individuals */ + /* wav and mw are no more changed */ + for(i=1; i<=imx; i++){ for(mi=1; mi
  • Result files (first order: no variance)\n \
  • Result files (second order (variance)\n \ "); - fprintf(fichtm,"

    • Result files (first order: no variance)

      \n \ - - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s
      \n ", - jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_")); + fprintf(fichtm,"

      • Result files (first order: no variance)

        \n"); + fprintf(fichtm,"
      • - Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s (html file)
        \n", + jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm")); + fprintf(fichtm,"
      • - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s (html file) ", + jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm")); + fprintf(fichtm,", %s (text file)
        \n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_")); fprintf(fichtm,"\ - Estimated transition probabilities over %d (stepm) months: %s
        \n ", stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_")); @@ -4822,12 +4990,14 @@ void printinghtml(char fileresu[], char - Period (stable) prevalence in each health state: %s
        \n", subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); fprintf(fichtm,"\ - - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ + - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ %s
        \n", estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); - fprintf(fichtm,"\ - - Population projections by age and states: \ + if(prevfcast==1){ + fprintf(fichtm,"\ + - Prevalence projections by age and states: \ %s
        \n
        • ", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); + } fprintf(fichtm," \n
        • Graphs

          • "); @@ -4847,16 +5017,16 @@ fprintf(fichtm," \n

          • Graphs fprintf(fichtm," ************\n
            "); } /* aij, bij */ - fprintf(fichtm,"
            - Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: %s_%d-1.svg
            \ -",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); + fprintf(fichtm,"
            - Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: %s_%d-1.svg
            \ +",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); /* Pij */ - fprintf(fichtm,"
            \n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s_%d-2.svg
            \ + fprintf(fichtm,"
            \n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s_%d-2.svg
            \ ",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); /* Quasi-incidences */ - fprintf(fichtm,"
            \n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ + fprintf(fichtm,"
            \n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too,\ - incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \ -divided by h: hPij/h : %s_%d-3.svg
            \ + incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \ +divided by h: hPij/h : %s_%d-3.svg
            \ ",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); /* Survival functions (period) in state j */ for(cpt=1; cpt<=nlstate;cpt++){ @@ -4874,6 +5044,14 @@ divided by h: hPij/h : \n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s%d_%d.svg
            \ ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); } + if(prevfcast==1){ + /* Projection of prevalence up to period (stable) prevalence in each health state */ + for(cpt=1; cpt<=nlstate;cpt++){ + fprintf(fichtm,"
            \n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s%d_%d.svg
            \ +", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1); + } + } + for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"\n
            - Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d%d.svg
            \ ",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); @@ -4961,10 +5139,11 @@ true period expectancies (those weighted } /******************* Gnuplot file **************/ -void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ + void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, char pathc[], double p[]){ char dirfileres[132],optfileres[132]; int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; + int lv=0, vlv=0, kl=0; int ng=0; int vpopbased; /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ @@ -4977,11 +5156,6 @@ void printinggnuplot(char fileresu[], ch /*#endif */ m=pow(2,cptcoveff); - /* Projected Prevalences */ -/* plot "NAGI0w_V1V2_monthlyb2b-proj/F_NAGI0w_V1V2_monthlyb2b-proj.txt" u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0))? $7/(1-$13):1/0) t 'p11' w line */ -/* replot "" u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0))? $8/(1-$14):1/0) t 'p21' w line */ -/* replot "" u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0)&&($9!=0))? $9/(1-$15):1/0) t 'p.1' w line */ - /* Contribution to likelihood */ /* Plot the probability implied in the likelihood */ fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); @@ -4993,15 +5167,15 @@ void printinggnuplot(char fileresu[], ch /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */ /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */ fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_")); - fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); + fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_")); - fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); + fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); for (i=1; i<= nlstate ; i ++) { fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i); fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk)); - fprintf(ficgp," u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1); + fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1); for (j=2; j<= nlstate+ndeath ; j ++) { - fprintf(ficgp,",\\\n \"\" u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j); + fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j); } fprintf(ficgp,";\nset out; unset ylabel;\n"); } @@ -5014,9 +5188,20 @@ void printinggnuplot(char fileresu[], ch strcpy(dirfileres,optionfilefiname); strcpy(optfileres,"vpl"); /* 1eme*/ - fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n"); - for (cpt=1; cpt<= nlstate ; cpt ++) { - for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ + for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ + for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */ + /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ + fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); + for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + fprintf(ficgp," V%d=%d ",k,vlv); + } + fprintf(ficgp,"\n#\n"); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1); fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \n\ @@ -5043,8 +5228,18 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u } /* k1 */ } /* cpt */ /*2 eme*/ - fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); for (k1=1; k1<= m ; k1 ++) { + fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); + for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + fprintf(ficgp," V%d=%d ",k,vlv); + } + fprintf(ficgp,"\n#\n"); + fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ if(vpopbased==0) @@ -5077,10 +5272,22 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u } /* vpopbased */ fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ } /* k1 */ + + /*3eme*/ - for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { + fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%d state=%d",k1, cpt); + for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + fprintf(ficgp," V%d=%d ",k,vlv); + } + fprintf(ficgp,"\n#\n"); + /* k=2+nlstate*(2*cpt-2); */ k=2+(nlstate+1)*(cpt-1); fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); @@ -5106,13 +5313,23 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u /* Survival functions (period) from state i in state j by initial state i */ for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ - k=3; - fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt); + for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + fprintf(ficgp," V%d=%d ",k,vlv); + } + fprintf(ficgp,"\n#\n"); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ set ter svg size 640, 480\n\ unset log y\n\ plot [%.f:%.f] ", ageminpar, agemaxpar); + k=3; for (i=1; i<= nlstate ; i ++){ if(i==1) fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); @@ -5131,13 +5348,23 @@ plot [%.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 if any */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ - k=3; 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); + for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + fprintf(ficgp," V%d=%d ",k,vlv); + } + fprintf(ficgp,"\n#\n"); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ set ter svg size 640, 480\n\ unset log y\n\ plot [%.f:%.f] ", ageminpar, agemaxpar); + k=3; for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ if(j==1) fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); @@ -5163,15 +5390,25 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) } /* end covariate */ /* CV preval stable (period) for each covariate */ - for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ + for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ - k=3; - fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); + for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + fprintf(ficgp," V%d=%d ",k,vlv); + } + fprintf(ficgp,"\n#\n"); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ set ter svg size 640, 480\n\ unset log y\n\ plot [%.f:%.f] ", ageminpar, agemaxpar); + k=3; /* Offset */ for (i=1; i<= nlstate ; i ++){ if(i==1) fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); @@ -5187,6 +5424,84 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) } /* end cpt state*/ } /* end covariate */ + if(prevfcast==1){ + /* Projection from cross-sectional to stable (period) for each covariate */ + + for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ + 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); + for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + fprintf(ficgp," V%d=%d ",k,vlv); + } + fprintf(ficgp,"\n#\n"); + + fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ +set ter svg size 640, 480\n\ +unset log y\n\ +plot [%.f:%.f] ", ageminpar, agemaxpar); + for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ + /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ + /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ + if(i==1){ + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); + }else{ + fprintf(ficgp,",\\\n '' "); + } + if(cptcoveff ==0){ /* No covariate */ + fprintf(ficgp," u 2:("); /* Age is in 2 */ + /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ + if(i==nlstate+1) + fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \ + 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,cpt ); + else + fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ + 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); + }else{ + fprintf(ficgp,"u 6:(("); /* Age is in 6 */ + /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ + kl=0; + for (k=1; k<=cptcoveff; k++){ /* For each covariate */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[lv]][lv]; + kl++; + /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ + /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ + /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ + /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ + if(k==cptcoveff) + if(i==nlstate+1) + fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \ + 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,cpt ); + else + fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \ + 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); + else{ + fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]); + kl++; + } + } /* end covariate */ + } /* end if covariate */ + } /* nlstate */ + fprintf(ficgp,"\nset out\n"); + } /* end cpt state*/ + } /* end covariate */ + } /* End if prevfcast */ + + /* proba elementaires */ fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); for(i=1,jk=1; i <=nlstate; i++){ @@ -5378,6 +5693,12 @@ void prevforecast(char fileres[], double char fileresf[FILENAMELENGTH]; agelim=AGESUP; + /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people + in each health status at the date of interview (if between dateprev1 and dateprev2). + We still use firstpass and lastpass as another selection. + */ + /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */ + /* firstpass, lastpass, stepm, weightopt, model); */ prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); strcpy(fileresf,"F_"); @@ -5386,8 +5707,8 @@ void prevforecast(char fileres[], double printf("Problem with forecast resultfile: %s\n", fileresf); fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf); } - printf("Computing forecasting: result on file '%s' \n", fileresf); - fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf); + printf("Computing forecasting: result on file '%s', please wait... \n", fileresf); + fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf); if (cptcoveff==0) ncodemax[cptcoveff]=1; @@ -5428,12 +5749,11 @@ void prevforecast(char fileres[], double for(cptcov=1, k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ k=k+1; - fprintf(ficresf,"\n#******"); + fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } - fprintf(ficresf,"******\n"); - fprintf(ficresf,"# Covariate valuofcovar yearproj age"); + fprintf(ficresf," yearproj age"); for(j=1; j<=nlstate+ndeath;j++){ for(i=1; i<=nlstate;i++) fprintf(ficresf," p%d%d",i,j); @@ -5483,6 +5803,9 @@ void prevforecast(char fileres[], double if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fclose(ficresf); + printf("End of Computing forecasting \n"); + fprintf(ficlog,"End of Computing forecasting\n"); + } /************** Forecasting *****not tested NB*************/ @@ -6322,12 +6645,12 @@ int calandcheckages(int imx, int maxwav, for (i=1; i<=imx; i++) { agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); for(m=firstpass; (m<= lastpass); m++){ - if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ + if(s[m][i] >0 || s[m][i]==-1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ /* What if s[m][i]=-1 */ if (s[m][i] >= nlstate+1) { if(agedc[i]>0){ if((int)moisdc[i]!=99 && (int)andc[i]!=9999){ agev[m][i]=agedc[i]; - /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ + /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ }else { if ((int)andc[i]!=9999){ nbwarn++; @@ -6337,7 +6660,7 @@ int calandcheckages(int imx, int maxwav, } } } /* agedc > 0 */ - } + } /* end if */ else if(s[m][i] !=9){ /* Standard case, age in fractional years but with the precision of a month */ agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); @@ -6353,17 +6676,23 @@ int calandcheckages(int imx, int maxwav, } /*agev[m][i]=anint[m][i]-annais[i];*/ /* agev[m][i] = age[i]+2*m;*/ - } + } /* en if 9*/ else { /* =9 */ + /* printf("Debug num[%d]=%ld s[%d][%d]=%d\n",i,num[i], m,i, s[m][i]); */ agev[m][i]=1; s[m][i]=-1; } } - else /*= 0 Unknown */ + else if(s[m][i]==0) /*= 0 Unknown */ agev[m][i]=1; - } - + else{ + printf("Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]); + fprintf(ficlog, "Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]); + agev[m][i]=0; + } + } /* End for lastpass */ } + for (i=1; i<=imx; i++) { for(m=firstpass; (m<=lastpass); m++){ if (s[m][i] > (nlstate+ndeath)) { @@ -7362,19 +7691,32 @@ Please run with mle=-1 to get a correct free_vector(annais,1,n); /* free_matrix(mint,1,maxwav,1,n); free_matrix(anint,1,maxwav,1,n);*/ - free_vector(moisdc,1,n); - free_vector(andc,1,n); + /* free_vector(moisdc,1,n); */ + /* free_vector(andc,1,n); */ /* */ wav=ivector(1,imx); - dh=imatrix(1,lastpass-firstpass+1,1,imx); - bh=imatrix(1,lastpass-firstpass+1,1,imx); - mw=imatrix(1,lastpass-firstpass+1,1,imx); + /* dh=imatrix(1,lastpass-firstpass+1,1,imx); */ + /* bh=imatrix(1,lastpass-firstpass+1,1,imx); */ + /* mw=imatrix(1,lastpass-firstpass+1,1,imx); */ + dh=imatrix(1,lastpass-firstpass+2,1,imx); /* We are adding a wave if status is unknown at last wave but death occurs after last wave.*/ + bh=imatrix(1,lastpass-firstpass+2,1,imx); + mw=imatrix(1,lastpass-firstpass+2,1,imx); /* Concatenates waves */ + /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i. + Death is a valid wave (if date is known). + mw[mi][i] is the number of (mi=1 to wav[i]) effective wave out of mi of individual i + dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] + and mw[mi+1][i]. dh depends on stepm. + */ + concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); /* */ + free_vector(moisdc,1,n); + free_vector(andc,1,n); + /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); @@ -7382,13 +7724,20 @@ Please run with mle=-1 to get a correct Ndum =ivector(-1,NCOVMAX); if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ - /* Nbcode gives the value of the lth modality of jth covariate, in + /* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ - /* 1 to ncodemax[j] is the maximum value of this jth covariate */ + /* 1 to ncodemax[j] which is the maximum value of this jth covariate */ /* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/ /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ + /* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j, + * codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded + * (currently 0 or 1) in the data. + * In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of + * corresponding modality (h,j). + */ + h=0; @@ -7398,8 +7747,9 @@ Please run with mle=-1 to get a correct m=pow(2,cptcoveff); /**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 - * For k=4 covariates, h goes from 1 to 2**k - * codtabm(h,k)= 1 & (h-1) >> (k-1) ; + * For k=4 covariates, h goes from 1 to m=2**k + * codtabm(h,k)= (1 & (h-1) >> (k-1)) + 1; + * #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 * h\k 1 2 3 4 *______________________________ * 1 i=1 1 i=1 1 i=1 1 i=1 1 @@ -7419,6 +7769,49 @@ Please run with mle=-1 to get a correct * 15 i=8 1 2 2 2 * 16 2 2 2 2 */ + /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */ + /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4 + * and the value of each covariate? + * V1=1, V2=1, V3=2, V4=1 ? + * h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok. + * h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st. + * In order to get the real value in the data, we use nbcode + * nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0 + * We are keeping this crazy system in order to be able (in the future?) + * to have more than 2 values (0 or 1) for a covariate. + * #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 + * h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1 + * bbbbbbbb + * 76543210 + * h-1 00000101 (6-1=5) + *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift + * & + * 1 00000001 (1) + * 00000001 = 1 & ((h-1) >> (k-1)) + * +1= 00000010 =2 + * + * h=14, k=3 => h'=h-1=13, k'=k-1=2 + * h' 1101 =2^3+2^2+0x2^1+2^0 + * >>k' 11 + * & 00000001 + * = 00000001 + * +1 = 00000010=2 = codtabm(14,3) + * Reverse h=6 and m=16? + * cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1. + * for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff) + * decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1 + * decodtabm(h,j,cptcoveff)= (h <= (1<> (j-1)) & 1) +1 : -1) + * V3=decodtabm(14,3,2**4)=2 + * h'=13 1101 =2^3+2^2+0x2^1+2^0 + *(h-1) >> (j-1) 0011 =13 >> 2 + * &1 000000001 + * = 000000001 + * +1= 000000010 =2 + * 2211 + * V1=1+1, V2=0+1, V3=1+1, V4=1+1 + * V3=2 + */ + /* /\* for(h=1; h <=100 ;h++){ *\/ */ /* /\* printf("h=%2d ", h); *\/ */ /* /\* for(k=1; k <=10; k++){ *\/ */ @@ -7497,7 +7890,7 @@ Title=%s
            Datafile=%s Firstpass=%d La optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(fichtm,"\n\n\nIMaCh %s\n IMaCh for Interpolated Markov Chain
            \nSponsored by Copyright (C) 2002-2015 INED-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015
            \ + fprintf(fichtm,"\n\n\nIMaCh %s\n IMaCh for Interpolated Markov Chain
            \nSponsored by Copyright (C) 2002-2015 INED-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015-2018
            \
            \n\ IMaCh-%s
            %s             \
            \n\ @@ -7527,7 +7920,8 @@ Title=%s
            Datafile=%s Firstpass=%d La /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ + firstpass, lastpass, stepm, weightopt, model); fprintf(fichtm,"\n"); fprintf(fichtm,"
            Total number of observations=%d
            \n\ @@ -7990,8 +8384,8 @@ Please run with mle=-1 to get a correct } fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); - fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); - fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); + fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl); + fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl); /* Other stuffs, more or less useful */ while((c=getc(ficpar))=='#' && c!= EOF){ @@ -8054,19 +8448,19 @@ Please run with mle=-1 to get a correct This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); }else - printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); + printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p); printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\ - model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ - jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); + model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \ + jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2); /*------------ free_vector -------------*/ /* chdir(path); */ - free_ivector(wav,1,imx); - free_imatrix(dh,1,lastpass-firstpass+1,1,imx); - free_imatrix(bh,1,lastpass-firstpass+1,1,imx); - free_imatrix(mw,1,lastpass-firstpass+1,1,imx); + /* free_ivector(wav,1,imx); */ /* Moved after last prevalence call */ + /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ + /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ + /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ free_lvector(num,1,n); free_vector(agedc,1,n); /*free_matrix(covar,0,NCOVMAX,1,n);*/ @@ -8126,6 +8520,11 @@ Please run with mle=-1 to get a correct /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); */ + free_ivector(wav,1,imx); + free_imatrix(dh,1,lastpass-firstpass+2,1,imx); + free_imatrix(bh,1,lastpass-firstpass+2,1,imx); + free_imatrix(mw,1,lastpass-firstpass+2,1,imx); + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);