|
|
| version 1.210, 2015/11/18 17:41:20 | version 1.216, 2015/12/18 17:32:11 |
|---|---|
| Line 1 | Line 1 |
| /* $Id$ | /* $Id$ |
| $State$ | $State$ |
| $Log$ | $Log$ |
| Revision 1.216 2015/12/18 17:32:11 brouard | |
| Summary: 0.98r4 Warning and status=-2 | |
| Version 0.98r4 is now: | |
| - displaying an error when status is -1, date of interview unknown and date of death known; | |
| - permitting a status -2 when the vital status is unknown at a known date of right truncation. | |
| Older changes concerning s=-2, dating from 2005 have been supersed. | |
| 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 | Revision 1.210 2015/11/18 17:41:20 brouard |
| Summary: Start working on projected prevalences | Summary: Start working on projected prevalences |
| Line 757 typedef struct { | Line 782 typedef struct { |
| #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ | #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
| #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ | #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
| #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 | #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| /*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/ | |
| #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 | |
| #define MAXN 20000 | #define MAXN 20000 |
| #define YEARM 12. /**< Number of months per year */ | #define YEARM 12. /**< Number of months per year */ |
| #define AGESUP 130 | #define AGESUP 130 |
| Line 814 double **matprod2(); /* test */ | Line 841 double **matprod2(); /* test */ |
| double **oldm, **newm, **savm; /* Working pointers to matrices */ | double **oldm, **newm, **savm; /* Working pointers to matrices */ |
| double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ | double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ |
| /*FILE *fic ; */ /* Used in readdata only */ | /*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; | FILE *ficlog, *ficrespow; |
| int globpr=0; /* Global variable for printing or not */ | int globpr=0; /* Global variable for printing or not */ |
| double fretone; /* Only one call to likelihood */ | double fretone; /* Only one call to likelihood */ |
| Line 1376 char *subdirf3(char fileres[], char *pre | Line 1403 char *subdirf3(char fileres[], char *pre |
| strcat(tmpout,fileres); | strcat(tmpout,fileres); |
| return tmpout; | 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[]) | char *asc_diff_time(long time_sec, char ascdiff[]) |
| { | { |
| long sec_left, days, hours, minutes; | long sec_left, days, hours, minutes; |
| Line 2213 double ***hpxij(double ***po, int nhstep | Line 2263 double ***hpxij(double ***po, int nhstep |
| double **out, cov[NCOVMAX+1]; | double **out, cov[NCOVMAX+1]; |
| double **newm; | double **newm; |
| double agexact; | double agexact; |
| double agebegin, ageend; | |
| /* Hstepm could be zero and should return the unit matrix */ | /* Hstepm could be zero and should return the unit matrix */ |
| for (i=1;i<=nlstate+ndeath;i++) | for (i=1;i<=nlstate+ndeath;i++) |
| Line 2226 double ***hpxij(double ***po, int nhstep | Line 2277 double ***hpxij(double ***po, int nhstep |
| newm=savm; | newm=savm; |
| /* Covariates have to be included here again */ | /* Covariates have to be included here again */ |
| cov[1]=1.; | 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; | cov[2]=agexact; |
| if(nagesqr==1) | if(nagesqr==1) |
| cov[3]= agexact*agexact; | cov[3]= agexact*agexact; |
| Line 2403 double func( double *x) | Line 2454 double func( double *x) |
| /* else */ | /* else */ |
| /* lli=log(out[s1][s2] - savm[s1][s2]); */ | /* lli=log(out[s1][s2] - savm[s1][s2]); */ |
| /* #endif */ | /* #endif */ |
| lli=log(out[s1][s2] - savm[s1][s2]); | lli=log(out[s1][s2] - savm[s1][s2]); |
| } else if (s2==-2) { | } else if ( s2==-1 ) { /* alive */ |
| for (j=1,survp=0. ; j<=nlstate; j++) | for (j=1,survp=0. ; j<=nlstate; j++) |
| survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; | survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
| /*survp += out[s1][j]; */ | /*survp += out[s1][j]; */ |
| lli= log(survp); | lli= log(survp); |
| } | } |
| else if (s2==-4) { | else if (s2==-4) { |
| for (j=3,survp=0. ; j<=nlstate; j++) | for (j=3,survp=0. ; j<=nlstate; j++) |
| survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; | survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
| lli= log(survp); | lli= log(survp); |
| } | } |
| else if (s2==-5) { | else if (s2==-5) { |
| for (j=1,survp=0. ; j<=2; j++) | for (j=1,survp=0. ; j<=2; j++) |
| survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; | survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
| lli= log(survp); | lli= log(survp); |
| } | } |
| else{ | else{ |
| lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ | lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
| /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ | /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ |
| Line 2535 double func( double *x) | Line 2583 double func( double *x) |
| s2=s[mw[mi+1][i]][i]; | s2=s[mw[mi+1][i]][i]; |
| if( s2 > nlstate){ | if( s2 > nlstate){ |
| lli=log(out[s1][s2] - savm[s1][s2]); | lli=log(out[s1][s2] - savm[s1][s2]); |
| } else if ( s2==-1 ) { /* alive */ | |
| for (j=1,survp=0. ; j<=nlstate; j++) | |
| survp += out[s1][j]; | |
| lli= log(survp); | |
| }else{ | }else{ |
| lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ | lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
| } | } |
| Line 2597 double funcone( double *x) | Line 2649 double funcone( double *x) |
| int s1, s2; | int s1, s2; |
| double bbh, survp; | double bbh, survp; |
| double agexact; | double agexact; |
| double agebegin, ageend; | |
| /*extern weight */ | /*extern weight */ |
| /* We are differentiating ll according to initial status */ | /* We are differentiating ll according to initial status */ |
| /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ | /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
| Line 2615 double funcone( double *x) | Line 2668 double funcone( double *x) |
| oldm[ii][j]=(ii==j ? 1.0 : 0.0); | oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
| savm[ii][j]=(ii==j ? 1.0 : 0.0); | savm[ii][j]=(ii==j ? 1.0 : 0.0); |
| } | } |
| for(d=0; d<dh[mi][i]; d++){ | |
| agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ | |
| ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ | |
| for(d=0; d<dh[mi][i]; d++){ /* Delay between two effective waves */ | |
| /*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.*/ | |
| newm=savm; | newm=savm; |
| agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; | agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
| cov[2]=agexact; | cov[2]=agexact; |
| Line 2636 double funcone( double *x) | Line 2694 double funcone( double *x) |
| s1=s[mw[mi][i]][i]; | s1=s[mw[mi][i]][i]; |
| s2=s[mw[mi+1][i]][i]; | s2=s[mw[mi+1][i]][i]; |
| if(s2==-1){ | |
| printf(" s1=%d, s2=%d i=%d \n", s1, s2, i); | |
| /* exit(1); */ | |
| } | |
| bbh=(double)bh[mi][i]/(double)stepm; | bbh=(double)bh[mi][i]/(double)stepm; |
| /* bias is positive if real duration | /* bias is positive if real duration |
| * is higher than the multiple of stepm and negative otherwise. | * is higher than the multiple of stepm and negative otherwise. |
| */ | */ |
| if( s2 > nlstate && (mle <5) ){ /* Jackson */ | if( s2 > nlstate && (mle <5) ){ /* Jackson */ |
| lli=log(out[s1][s2] - savm[s1][s2]); | lli=log(out[s1][s2] - savm[s1][s2]); |
| } else if (s2==-2) { | } else if ( s2==-1 ) { /* alive */ |
| for (j=1,survp=0. ; j<=nlstate; j++) | for (j=1,survp=0. ; j<=nlstate; j++) |
| survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; | survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
| lli= log(survp); | lli= log(survp); |
| Line 2663 double funcone( double *x) | Line 2725 double funcone( double *x) |
| ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; | ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ | /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ |
| if(globpr){ | if(globpr){ |
| fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ | fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ |
| %11.6f %11.6f %11.6f ", \ | %11.6f %11.6f %11.6f ", \ |
| num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, | num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, |
| 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); | 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); |
| for(k=1,llt=0.,l=0.; k<=nlstate; k++){ | for(k=1,llt=0.,l=0.; k<=nlstate; k++){ |
| llt +=ll[k]*gipmx/gsw; | llt +=ll[k]*gipmx/gsw; |
| Line 2703 void likelione(FILE *ficres,double p[], | Line 2765 void likelione(FILE *ficres,double p[], |
| printf("Problem with resultfile: %s\n", fileresilk); | printf("Problem with resultfile: %s\n", fileresilk); |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); | fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); |
| } | } |
| fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); | fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); |
| fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav "); | fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav "); |
| /* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ | /* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ |
| for(k=1; k<=nlstate; k++) | for(k=1; k<=nlstate; k++) |
| fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); | fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); |
| Line 2722 void likelione(FILE *ficres,double p[], | Line 2784 void likelione(FILE *ficres,double p[], |
| for (k=1; k<= nlstate ; k++) { | for (k=1; k<= nlstate ; k++) { |
| fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ | fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ |
| <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); | <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); |
| } | } |
| fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \ | fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \ |
| Line 3178 void pstamp(FILE *fichier) | Line 3240 void pstamp(FILE *fichier) |
| } | } |
| /************ Frequencies ********************/ | /************ 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 */ | { /* Some frequencies */ |
| int i, m, jk, j1, bool, z1,j; | int i, m, jk, j1, bool, z1,j; |
| int mi; /* Effective wave */ | |
| int first; | int first; |
| double ***freq; /* Frequencies */ | double ***freq; /* Frequencies */ |
| double *pp, **prop; | double *pp, **prop; |
| double pos,posprop, k2, dateintsum=0,k2cpt=0; | 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); | pp=vector(1,nlstate); |
| prop=matrix(1,nlstate,iagemin,iagemax+3); | prop=matrix(1,nlstate,iagemin,iagemax+3); |
| strcpy(fileresp,"P_"); | strcpy(fileresp,"P_"); |
| strcat(fileresp,fileresu); | strcat(fileresp,fileresu); |
| /*strcat(fileresphtm,fileresu);*/ | |
| if((ficresp=fopen(fileresp,"w"))==NULL) { | if((ficresp=fopen(fileresp,"w"))==NULL) { |
| printf("Problem with prevalence resultfile: %s\n", fileresp); | printf("Problem with prevalence resultfile: %s\n", fileresp); |
| fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); | fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); |
| exit(0); | 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,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ | |
| <hr size=\"2\" color=\"#EC5E5E\"> \n\ | |
| Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ | |
| fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); | |
| } | |
| fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\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,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ | |
| <hr size=\"2\" color=\"#EC5E5E\"> \n\ | |
| Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ | |
| fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); | |
| } | |
| fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr); | |
| freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); | freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); |
| j1=0; | j1=0; |
| Line 3205 void freqsummary(char fileres[], int ia | Line 3303 void freqsummary(char fileres[], int ia |
| first=1; | first=1; |
| /* for(k1=1; k1<=j ; k1++){ */ /* Loop on covariates */ | for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */ |
| /* for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */ | |
| /* j1++; */ | |
| for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ | |
| /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); | /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
| scanf("%d", i);*/ | scanf("%d", i);*/ |
| for (i=-5; i<=nlstate+ndeath; i++) | for (i=-5; i<=nlstate+ndeath; i++) |
| Line 3222 void freqsummary(char fileres[], int ia | Line 3317 void freqsummary(char fileres[], int ia |
| dateintsum=0; | dateintsum=0; |
| k2cpt=0; | k2cpt=0; |
| for (i=1; i<=imx; i++) { | for (i=1; i<=imx; i++) { /* For each individual i */ |
| bool=1; | bool=1; |
| if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ | if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
| for (z1=1; z1<=cptcoveff; z1++) | for (z1=1; z1<=cptcoveff; z1++) |
| Line 3235 void freqsummary(char fileres[], int ia | Line 3330 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*/ | /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
| } | } |
| } /* cptcovn > 0 */ | } /* cptcovn > 0 */ |
| if (bool==1){ | if (bool==1){ |
| for(m=firstpass; m<=lastpass; m++){ | /* for(m=firstpass; m<=lastpass; m++){ */ |
| k2=anint[m][i]+(mint[m][i]/12.); | for(mi=1; mi<wav[i];mi++){ |
| /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ | m=mw[mi][i]; |
| if(agev[m][i]==0) agev[m][i]=iagemax+1; | /* dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective (mi) waves m=mw[mi][i] |
| if(agev[m][i]==1) agev[m][i]=iagemax+2; | and mw[mi+1][i]. dh depends on stepm. */ |
| if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; | agebegin=agev[m][i]; /* Age at beginning of wave before transition*/ |
| ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /* Age at end of wave and transition */ | |
| if(m >=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 (m<lastpass) { | if (m<lastpass) { |
| freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; | /* if(s[m][i]==4 && s[m+1][i]==4) */ |
| freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; | /* 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)) { | if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) { |
| dateintsum=dateintsum+k2; | dateintsum=dateintsum+k2; |
| k2cpt++; | k2cpt++; |
| /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ | /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ |
| } | } |
| /*}*/ | /*}*/ |
| } /* end m */ | } /* end m */ |
| } /* end bool */ | } /* end bool */ |
| } /* end i = 1 to imx */ | } /* end i = 1 to imx */ |
| Line 3262 void freqsummary(char fileres[], int ia | Line 3370 void freqsummary(char fileres[], int ia |
| pstamp(ficresp); | pstamp(ficresp); |
| if (cptcovn>0) { | if (cptcovn>0) { |
| fprintf(ficresp, "\n#********** Variable "); | fprintf(ficresp, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); | fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
| fprintf(ficresp, "**********\n#"); | fprintf(ficresphtmfr, "\n<br/><br/><h3>********** 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, "**********</h3>\n"); | |
| fprintf(ficresphtmfr, "**********</h3>\n"); | |
| fprintf(ficlog, "\n#********** Variable "); | fprintf(ficlog, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); | 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,"<table style=\"text-align:center; border: 1px solid\">"); |
| for(i=1; i<=nlstate;i++) { | |
| fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); | fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); |
| fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); | |
| } | |
| fprintf(ficresp, "\n"); | fprintf(ficresp, "\n"); |
| fprintf(ficresphtm, "\n"); | |
| /* Header of frequency table by age */ | |
| fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">"); | |
| fprintf(ficresphtmfr,"<th>Age</th> "); | |
| for(jk=-1; jk <=nlstate+ndeath; jk++){ | |
| for(m=-1; m <=nlstate+ndeath; m++){ | |
| if(jk!=0 && m!=0) | |
| fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m); | |
| } | |
| } | |
| fprintf(ficresphtmfr, "\n"); | |
| /* For each age */ | |
| for(i=iagemin; i <= iagemax+3; i++){ | for(i=iagemin; i <= iagemax+3; i++){ |
| if(i==iagemax+3){ | fprintf(ficresphtm,"<tr>"); |
| if(i==iagemax+1){ | |
| fprintf(ficlog,"1"); | |
| fprintf(ficresphtmfr,"<tr><th>0</th> "); | |
| }else if(i==iagemax+2){ | |
| fprintf(ficlog,"0"); | |
| fprintf(ficresphtmfr,"<tr><th>Unknown</th> "); | |
| }else if(i==iagemax+3){ | |
| fprintf(ficlog,"Total"); | fprintf(ficlog,"Total"); |
| fprintf(ficresphtmfr,"<tr><th>Total</th> "); | |
| }else{ | }else{ |
| if(first==1){ | if(first==1){ |
| first=0; | first=0; |
| printf("See log file for details...\n"); | printf("See log file for details...\n"); |
| } | } |
| fprintf(ficresphtmfr,"<tr><th>%d</th> ",i); | |
| fprintf(ficlog,"Age %d", i); | fprintf(ficlog,"Age %d", i); |
| } | } |
| for(jk=1; jk <=nlstate ; jk++){ | for(jk=1; jk <=nlstate ; jk++){ |
| Line 3322 void freqsummary(char fileres[], int ia | Line 3463 void freqsummary(char fileres[], int ia |
| if( i <= iagemax){ | if( i <= iagemax){ |
| if(pos>=1.e-5){ | if(pos>=1.e-5){ |
| fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); | fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); |
| fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",i,prop[jk][i]/posprop, prop[jk][i],posprop); | |
| /*probs[i][jk][j1]= pp[jk]/pos;*/ | /*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]);*/ | /*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(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); |
| fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",i, prop[jk][i],posprop); | |
| } | |
| } | } |
| } | } |
| for(jk=-1; jk <=nlstate+ndeath; jk++) | for(jk=-1; jk <=nlstate+ndeath; jk++){ |
| for(m=-1; m <=nlstate+ndeath; m++) | for(m=-1; m <=nlstate+ndeath; m++){ |
| if(freq[jk][m][i] !=0 ) { | if(freq[jk][m][i] !=0 ) { /* minimizing output */ |
| if(first==1) | if(first==1){ |
| printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); | printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); |
| } | |
| fprintf(ficlog," %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,"<td>%.0f</td> ",freq[jk][m][i]); | |
| } | |
| } | |
| fprintf(ficresphtmfr,"</tr>\n "); | |
| if(i <= iagemax){ | |
| fprintf(ficresp,"\n"); | fprintf(ficresp,"\n"); |
| fprintf(ficresphtm,"</tr>\n"); | |
| } | |
| if(first==1) | if(first==1) |
| printf("Others in log...\n"); | printf("Others in log...\n"); |
| fprintf(ficlog,"\n"); | fprintf(ficlog,"\n"); |
| } /* end loop i */ | } /* end loop i */ |
| fprintf(ficresphtm,"</table>\n"); | |
| fprintf(ficresphtmfr,"</table>\n"); | |
| /*}*/ | /*}*/ |
| } /* end j1 */ | } /* end j1 */ |
| dateintmean=dateintsum/k2cpt; | dateintmean=dateintsum/k2cpt; |
| fclose(ficresp); | fclose(ficresp); |
| fclose(ficresphtm); | |
| fclose(ficresphtmfr); | |
| free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); | free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); |
| free_vector(pp,1,nlstate); | free_vector(pp,1,nlstate); |
| free_matrix(prop,1,nlstate,iagemin, iagemax+3); | free_matrix(prop,1,nlstate,iagemin, iagemax+3); |
| Line 3363 void prevalence(double ***probs, double | Line 3519 void prevalence(double ***probs, double |
| */ | */ |
| int i, m, jk, j1, bool, z1,j; | int i, m, jk, j1, bool, z1,j; |
| int mi; /* Effective wave */ | |
| int iage; | |
| double agebegin, ageend; | |
| double **prop; | double **prop; |
| double posprop; | double posprop; |
| Line 3382 void prevalence(double ***probs, double | Line 3541 void prevalence(double ***probs, double |
| first=1; | first=1; |
| for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ | for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ |
| /*for(i1=1; i1<=ncodemax[k1];i1++){ | for (i=1; i<=nlstate; i++) |
| j1++;*/ | for(iage=iagemin; iage <= iagemax+3; iage++) |
| prop[i][iage]=0.0; | |
| for (i=1; i<=nlstate; i++) | |
| for(m=iagemin; m <= iagemax+3; m++) | for (i=1; i<=imx; i++) { /* Each individual */ |
| prop[i][m]=0.0; | bool=1; |
| if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ | |
| for (i=1; i<=imx; i++) { /* Each individual */ | for (z1=1; z1<=cptcoveff; z1++) |
| bool=1; | if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) |
| if (cptcovn>0) { | bool=0; |
| for (z1=1; z1<=cptcoveff; z1++) | } |
| if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) | if (bool==1) { |
| bool=0; | /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ |
| } | for(mi=1; mi<wav[i];mi++){ |
| if (bool==1) { | m=mw[mi][i]; |
| for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/ | agebegin=agev[m][i]; /* Age at beginning of wave before transition*/ |
| /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */ | |
| if(m >=firstpass && m <=lastpass){ | |
| y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ | 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 ((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]==0) agev[m][i]=iagemax+1; |
| if(agev[m][i]==1) agev[m][i]=iagemax+2; | if(agev[m][i]==1) agev[m][i]=iagemax+2; |
| if((int)agev[m][i] <iagemin || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); | if((int)agev[m][i] <iagemin || (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]]);*/ | /*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]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ |
| prop[s[m][i]][iagemax+3] += weight[i]; | prop[s[m][i]][iagemax+3] += weight[i]; |
| } | } /* end valid statuses */ |
| } | } /* end selection of dates */ |
| } /* end selection of waves */ | } /* end selection of waves */ |
| } | } /* end effective waves */ |
| } | } /* end bool */ |
| for(i=iagemin; i <= iagemax+3; i++){ | } |
| for(jk=1,posprop=0; jk <=nlstate ; jk++) { | for(i=iagemin; i <= iagemax+3; i++){ |
| posprop += prop[jk][i]; | for(jk=1,posprop=0; jk <=nlstate ; jk++) { |
| } | posprop += prop[jk][i]; |
| } | |
| for(jk=1; jk <=nlstate ; jk++){ | |
| if( i <= iagemax){ | for(jk=1; jk <=nlstate ; jk++){ |
| if(posprop>=1.e-5){ | if( i <= iagemax){ |
| probs[i][jk][j1]= prop[jk][i]/posprop; | if(posprop>=1.e-5){ |
| } else{ | probs[i][jk][j1]= prop[jk][i]/posprop; |
| if(first==1){ | } else{ |
| first=0; | if(first==1){ |
| 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]); | 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 i1 */ |
| } /* end j1 */ | } /* end j1 */ |
| Line 3452 void concatwav(int wav[], int **dh, int | Line 3614 void concatwav(int wav[], int **dh, int |
| int i, mi, m; | int i, mi, m; |
| /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; | /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; |
| double sum=0., jmean=0.;*/ | double sum=0., jmean=0.;*/ |
| int first; | int first, firstwo; |
| int j, k=0,jk, ju, jl; | int j, k=0,jk, ju, jl; |
| double sum=0.; | double sum=0.; |
| first=0; | first=0; |
| firstwo=0; | |
| jmin=100000; | jmin=100000; |
| jmax=-1; | jmax=-1; |
| jmean=0.; | jmean=0.; |
| for(i=1; i<=imx; i++){ | for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ |
| mi=0; | mi=0; |
| m=firstpass; | 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) | if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ |
| mw[++mi][i]=m; | mw[++mi][i]=m; |
| if(m >=lastpass) | } |
| if(m >=lastpass){ | |
| if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ | |
| printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); | |
| fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); | |
| mw[++mi][i]=m; | |
| } | |
| if(s[m][i]==-2){ /* Vital status is really unknown */ | |
| nbwarn++; | |
| if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ | |
| printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); | |
| fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); | |
| } | |
| break; | |
| } | |
| break; | break; |
| } | |
| else | else |
| m++; | m++; |
| }/* end while */ | }/* end while */ |
| if (s[m][i] > nlstate){ | |
| /* After last pass */ | |
| if (s[m][i] > nlstate){ /* In a death state */ | |
| mi++; /* Death is another wave */ | mi++; /* Death is another wave */ |
| /* if(mi==0) never been interviewed correctly before death */ | /* if(mi==0) never been interviewed correctly before death */ |
| /* Only death is a correct wave */ | /* Only death is a correct wave */ |
| mw[mi][i]=m; | mw[mi][i]=m; |
| }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */ | |
| /* m++; */ | |
| /* mi++; */ | |
| /* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */ | |
| /* mw[mi][i]=m; */ | |
| nberr++; | |
| if(firstwo==0){ | |
| printf("Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); | |
| fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); | |
| firstwo=1; | |
| }else if(firstwo==1){ | |
| fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); | |
| } | |
| } | } |
| wav[i]=mi; | wav[i]=mi; |
| if(mi==0){ | if(mi==0){ |
| nbwarn++; | nbwarn++; |
| Line 3489 void concatwav(int wav[], int **dh, int | Line 3681 void concatwav(int wav[], int **dh, int |
| } | } |
| } /* end mi==0 */ | } /* end mi==0 */ |
| } /* End individuals */ | } /* End individuals */ |
| /* wav and mw are no more changed */ | |
| for(i=1; i<=imx; i++){ | for(i=1; i<=imx; i++){ |
| for(mi=1; mi<wav[i];mi++){ | for(mi=1; mi<wav[i];mi++){ |
| if (stepm <=0) | if (stepm <=0) |
| Line 4804 To be simple, these graphs help to under | Line 4998 To be simple, these graphs help to under |
| void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ | void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ |
| int lastpass, int stepm, int weightopt, char model[],\ | int lastpass, int stepm, int weightopt, char model[],\ |
| int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ | int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ |
| int popforecast, int estepm ,\ | int popforecast, int prevfcast, int estepm , \ |
| double jprev1, double mprev1,double anprev1, \ | double jprev1, double mprev1,double anprev1, double dateprev1, \ |
| double jprev2, double mprev2,double anprev2){ | double jprev2, double mprev2,double anprev2, double dateprev2){ |
| int jj1, k1, i1, cpt; | int jj1, k1, i1, cpt; |
| fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ | fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
| <li><a href='#secondorder'>Result files (second order (variance)</a>\n \ | <li><a href='#secondorder'>Result files (second order (variance)</a>\n \ |
| </ul>"); | </ul>"); |
| fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \ | fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n"); |
| - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ", | fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n", |
| jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_")); | jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm")); |
| fprintf(fichtm,"<li> - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file) ", | |
| jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm")); | |
| fprintf(fichtm,", <a href=\"%s\">%s</a> (text file) <br>\n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_")); | |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", | - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", |
| stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_")); | stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_")); |
| Line 4822 void printinghtml(char fileresu[], char | Line 5019 void printinghtml(char fileresu[], char |
| - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", | - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); | subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); |
| fprintf(fichtm,"\ | 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, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . 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 href=\"%s\">%s</a> <br>\n", | <a href=\"%s\">%s</a> <br>\n", |
| estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); | estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); |
| fprintf(fichtm,"\ | if(prevfcast==1){ |
| - Population projections by age and states: \ | fprintf(fichtm,"\ |
| - Prevalence projections by age and states: \ | |
| <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); | <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); |
| } | |
| fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); | fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
| Line 4847 fprintf(fichtm," \n<ul><li><b>Graphs</b> | Line 5046 fprintf(fichtm," \n<ul><li><b>Graphs</b> |
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); | fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| } | } |
| /* aij, bij */ | /* aij, bij */ |
| fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \ | fprintf(fichtm,"<br>- 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: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \ |
| <img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); | <img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Pij */ | /* Pij */ |
| fprintf(fichtm,"<br>\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \ | fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \ |
| <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); | <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Quasi-incidences */ | /* Quasi-incidences */ |
| fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ | fprintf(fichtm,"<br>\n- I<sub>ij</sub> 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,\ | 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 \ | incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \ |
| divided by h: hPij/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \ | divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \ |
| <img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); | <img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Survival functions (period) in state j */ | /* Survival functions (period) in state j */ |
| for(cpt=1; cpt<=nlstate;cpt++){ | for(cpt=1; cpt<=nlstate;cpt++){ |
| Line 4874 divided by h: hPij/h : <a href=\"%s_%d-3 | Line 5073 divided by h: hPij/h : <a href=\"%s_%d-3 |
| fprintf(fichtm,"<br>\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. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ | fprintf(fichtm,"<br>\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. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ |
| <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); | <img src=\"%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,"<br>\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. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ | |
| <img src=\"%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++) { | for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"\n<br>- 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): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \ | fprintf(fichtm,"\n<br>- 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): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \ |
| <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); | <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); |
| Line 4961 true period expectancies (those weighted | Line 5168 true period expectancies (those weighted |
| } | } |
| /******************* Gnuplot file **************/ | /******************* 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]; | 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 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 ng=0; |
| int vpopbased; | int vpopbased; |
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ | /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
| Line 4977 void printinggnuplot(char fileresu[], ch | Line 5185 void printinggnuplot(char fileresu[], ch |
| /*#endif */ | /*#endif */ |
| m=pow(2,cptcoveff); | 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 */ | /* Contribution to likelihood */ |
| /* Plot the probability implied in the 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"); | fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); |
| Line 4993 void printinggnuplot(char fileresu[], ch | Line 5196 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)" */ | /* 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.svg\";",subdirf2(optionfilefiname,"ILK_")); */ |
| fprintf(ficgp,"\nset out \"%s-dest.png\";",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 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 ++) { | 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,"\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,"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 ++) { | 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"); | fprintf(ficgp,";\nset out; unset ylabel;\n"); |
| } | } |
| Line 5014 void printinggnuplot(char fileresu[], ch | Line 5217 void printinggnuplot(char fileresu[], ch |
| strcpy(dirfileres,optionfilefiname); | strcpy(dirfileres,optionfilefiname); |
| strcpy(optfileres,"vpl"); | strcpy(optfileres,"vpl"); |
| /* 1eme*/ | /* 1eme*/ |
| fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n"); | for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ |
| for (cpt=1; cpt<= nlstate ; cpt ++) { | for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */ |
| for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ | /* 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,"\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,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \n\ | fprintf(ficgp,"set xlabel \"Age\" \n\ |
| Line 5043 plot [%.f:%.f] \"%s\" every :::%d::%d u | Line 5257 plot [%.f:%.f] \"%s\" every :::%d::%d u |
| } /* k1 */ | } /* k1 */ |
| } /* cpt */ | } /* cpt */ |
| /*2 eme*/ | /*2 eme*/ |
| fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); | |
| for (k1=1; k1<= m ; k1 ++) { | 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); | 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*/ | for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| if(vpopbased==0) | if(vpopbased==0) |
| Line 5077 plot [%.f:%.f] \"%s\" every :::%d::%d u | Line 5301 plot [%.f:%.f] \"%s\" every :::%d::%d u |
| } /* vpopbased */ | } /* vpopbased */ |
| fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ | fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ |
| } /* k1 */ | } /* k1 */ |
| /*3eme*/ | /*3eme*/ |
| for (k1=1; k1<= m ; k1 ++) { | for (k1=1; k1<= m ; k1 ++) { |
| for (cpt=1; cpt<= nlstate ; cpt ++) { | 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*(2*cpt-2); */ |
| k=2+(nlstate+1)*(cpt-1); | k=2+(nlstate+1)*(cpt-1); |
| fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); | fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); |
| Line 5106 plot [%.f:%.f] \"%s\" every :::%d::%d u | Line 5342 plot [%.f:%.f] \"%s\" every :::%d::%d u |
| /* Survival functions (period) from state i in state j by initial state i */ | /* 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 (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ | 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,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ | fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
| set ter svg size 640, 480\n\ | set ter svg size 640, 480\n\ |
| unset log y\n\ | unset log y\n\ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); | plot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; | |
| for (i=1; i<= nlstate ; i ++){ | for (i=1; i<= nlstate ; i ++){ |
| if(i==1) | if(i==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); | fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
| Line 5131 plot [%.f:%.f] ", ageminpar, agemaxpar) | Line 5377 plot [%.f:%.f] ", ageminpar, agemaxpar) |
| /* Survival functions (period) from state i in state j by final state j */ | /* 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 (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ | 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); | 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,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ | fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
| set ter svg size 640, 480\n\ | set ter svg size 640, 480\n\ |
| unset log y\n\ | unset log y\n\ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); | plot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; | |
| for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ | for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ |
| if(j==1) | if(j==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); | fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
| Line 5163 plot [%.f:%.f] ", ageminpar, agemaxpar) | Line 5419 plot [%.f:%.f] ", ageminpar, agemaxpar) |
| } /* end covariate */ | } /* end covariate */ |
| /* CV preval stable (period) for each 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 */ | for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| k=3; | fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); |
| fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' 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,"P_"),cpt,k1); | fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ | fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
| set ter svg size 640, 480\n\ | set ter svg size 640, 480\n\ |
| unset log y\n\ | unset log y\n\ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); | plot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; /* Offset */ | |
| for (i=1; i<= nlstate ; i ++){ | for (i=1; i<= nlstate ; i ++){ |
| if(i==1) | if(i==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); | fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
| Line 5187 plot [%.f:%.f] ", ageminpar, agemaxpar) | Line 5453 plot [%.f:%.f] ", ageminpar, agemaxpar) |
| } /* end cpt state*/ | } /* end cpt state*/ |
| } /* end covariate */ | } /* 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 */ | /* proba elementaires */ |
| fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); | fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); |
| for(i=1,jk=1; i <=nlstate; i++){ | for(i=1,jk=1; i <=nlstate; i++){ |
| Line 5378 void prevforecast(char fileres[], double | Line 5722 void prevforecast(char fileres[], double |
| char fileresf[FILENAMELENGTH]; | char fileresf[FILENAMELENGTH]; |
| agelim=AGESUP; | 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); | prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
| strcpy(fileresf,"F_"); | strcpy(fileresf,"F_"); |
| Line 5386 void prevforecast(char fileres[], double | Line 5736 void prevforecast(char fileres[], double |
| printf("Problem with forecast resultfile: %s\n", fileresf); | printf("Problem with forecast resultfile: %s\n", fileresf); |
| fprintf(ficlog,"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); | printf("Computing forecasting: result on file '%s', please wait... \n", fileresf); |
| fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf); | fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf); |
| if (cptcoveff==0) ncodemax[cptcoveff]=1; | if (cptcoveff==0) ncodemax[cptcoveff]=1; |
| Line 5428 void prevforecast(char fileres[], double | Line 5778 void prevforecast(char fileres[], double |
| for(cptcov=1, k=0;cptcov<=i1;cptcov++){ | for(cptcov=1, k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ | for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ |
| k=k+1; | 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++) { | 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," yearproj age"); |
| fprintf(ficresf,"# Covariate valuofcovar yearproj age"); | |
| for(j=1; j<=nlstate+ndeath;j++){ | for(j=1; j<=nlstate+ndeath;j++){ |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| fprintf(ficresf," p%d%d",i,j); | fprintf(ficresf," p%d%d",i,j); |
| Line 5483 void prevforecast(char fileres[], double | Line 5832 void prevforecast(char fileres[], double |
| if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); | if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| fclose(ficresf); | fclose(ficresf); |
| printf("End of Computing forecasting \n"); | |
| fprintf(ficlog,"End of Computing forecasting\n"); | |
| } | } |
| /************** Forecasting *****not tested NB*************/ | /************** Forecasting *****not tested NB*************/ |
| Line 6302 int calandcheckages(int imx, int maxwav, | Line 6654 int calandcheckages(int imx, int maxwav, |
| for(m=2; (m<= maxwav); m++) { | for(m=2; (m<= maxwav); m++) { |
| if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){ | if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){ |
| anint[m][i]=9999; | anint[m][i]=9999; |
| s[m][i]=-1; | if (s[m][i] != -2) /* Keeping initial status of unknown vital status */ |
| s[m][i]=-1; | |
| } | } |
| if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){ | if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){ |
| *nberr = *nberr + 1; | *nberr = *nberr + 1; |
| Line 6322 int calandcheckages(int imx, int maxwav, | Line 6675 int calandcheckages(int imx, int maxwav, |
| for (i=1; i<=imx; i++) { | for (i=1; i<=imx; i++) { |
| agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); | agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); |
| for(m=firstpass; (m<= lastpass); m++){ | 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 (s[m][i] >= nlstate+1) { |
| if(agedc[i]>0){ | if(agedc[i]>0){ |
| if((int)moisdc[i]!=99 && (int)andc[i]!=9999){ | if((int)moisdc[i]!=99 && (int)andc[i]!=9999){ |
| agev[m][i]=agedc[i]; | 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 { | }else { |
| if ((int)andc[i]!=9999){ | if ((int)andc[i]!=9999){ |
| nbwarn++; | nbwarn++; |
| Line 6337 int calandcheckages(int imx, int maxwav, | Line 6690 int calandcheckages(int imx, int maxwav, |
| } | } |
| } | } |
| } /* agedc > 0 */ | } /* agedc > 0 */ |
| } | } /* end if */ |
| else if(s[m][i] !=9){ /* Standard case, age in fractional | else if(s[m][i] !=9){ /* Standard case, age in fractional |
| years but with the precision of a month */ | 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]); | agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); |
| Line 6353 int calandcheckages(int imx, int maxwav, | Line 6706 int calandcheckages(int imx, int maxwav, |
| } | } |
| /*agev[m][i]=anint[m][i]-annais[i];*/ | /*agev[m][i]=anint[m][i]-annais[i];*/ |
| /* agev[m][i] = age[i]+2*m;*/ | /* agev[m][i] = age[i]+2*m;*/ |
| } | } /* en if 9*/ |
| else { /* =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; | agev[m][i]=1; |
| s[m][i]=-1; | s[m][i]=-1; |
| } | } |
| } | } |
| else /*= 0 Unknown */ | else if(s[m][i]==0) /*= 0 Unknown */ |
| agev[m][i]=1; | 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 (i=1; i<=imx; i++) { |
| for(m=firstpass; (m<=lastpass); m++){ | for(m=firstpass; (m<=lastpass); m++){ |
| if (s[m][i] > (nlstate+ndeath)) { | if (s[m][i] > (nlstate+ndeath)) { |
| Line 7362 Please run with mle=-1 to get a correct | Line 7721 Please run with mle=-1 to get a correct |
| free_vector(annais,1,n); | free_vector(annais,1,n); |
| /* free_matrix(mint,1,maxwav,1,n); | /* free_matrix(mint,1,maxwav,1,n); |
| free_matrix(anint,1,maxwav,1,n);*/ | free_matrix(anint,1,maxwav,1,n);*/ |
| free_vector(moisdc,1,n); | /* free_vector(moisdc,1,n); */ |
| free_vector(andc,1,n); | /* free_vector(andc,1,n); */ |
| /* */ | /* */ |
| wav=ivector(1,imx); | wav=ivector(1,imx); |
| dh=imatrix(1,lastpass-firstpass+1,1,imx); | /* dh=imatrix(1,lastpass-firstpass+1,1,imx); */ |
| bh=imatrix(1,lastpass-firstpass+1,1,imx); | /* bh=imatrix(1,lastpass-firstpass+1,1,imx); */ |
| mw=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 */ |
| /* 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); | 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 */ | /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ |
| nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); | nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
| Line 7382 Please run with mle=-1 to get a correct | Line 7754 Please run with mle=-1 to get a correct |
| Ndum =ivector(-1,NCOVMAX); | Ndum =ivector(-1,NCOVMAX); |
| if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ | 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]; */ | 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).*/ | 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) */ | /* 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));*/ | /*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).*/ | /* 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; | h=0; |
| Line 7398 Please run with mle=-1 to get a correct | Line 7777 Please run with mle=-1 to get a correct |
| m=pow(2,cptcoveff); | m=pow(2,cptcoveff); |
| /**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 | /**< 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 | * For k=4 covariates, h goes from 1 to m=2**k |
| * codtabm(h,k)= 1 & (h-1) >> (k-1) ; | * 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 | * h\k 1 2 3 4 |
| *______________________________ | *______________________________ |
| * 1 i=1 1 i=1 1 i=1 1 i=1 1 | * 1 i=1 1 i=1 1 i=1 1 i=1 1 |
| Line 7419 Please run with mle=-1 to get a correct | Line 7799 Please run with mle=-1 to get a correct |
| * 15 i=8 1 2 2 2 | * 15 i=8 1 2 2 2 |
| * 16 2 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<<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++){ *\/ */ | /* /\* for(h=1; h <=100 ;h++){ *\/ */ |
| /* /\* printf("h=%2d ", h); *\/ */ | /* /\* printf("h=%2d ", h); *\/ */ |
| /* /\* for(k=1; k <=10; k++){ *\/ */ | /* /\* for(k=1; k <=10; k++){ *\/ */ |
| Line 7497 Title=%s <br>Datafile=%s Firstpass=%d La | Line 7920 Title=%s <br>Datafile=%s Firstpass=%d La |
| optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); | optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
| } | } |
| fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015</a></font><br> \ | fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \ |
| <hr size=\"2\" color=\"#EC5E5E\"> \n\ | <hr size=\"2\" color=\"#EC5E5E\"> \n\ |
| <font size=\"2\">IMaCh-%s <br> %s</font> \ | <font size=\"2\">IMaCh-%s <br> %s</font> \ |
| <hr size=\"2\" color=\"#EC5E5E\"> \n\ | <hr size=\"2\" color=\"#EC5E5E\"> \n\ |
| Line 7527 Title=%s <br>Datafile=%s Firstpass=%d La | Line 7950 Title=%s <br>Datafile=%s Firstpass=%d La |
| /* Calculates basic frequencies. Computes observed prevalence at single age | /* Calculates basic frequencies. Computes observed prevalence at single age |
| and prints on file fileres'p'. */ | 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,"\n"); |
| fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ | fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
| Line 7990 Please run with mle=-1 to get a correct | Line 8414 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 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(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\n",ageminpar,agemaxpar,bage,fage, estepm); | 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 */ | /* Other stuffs, more or less useful */ |
| while((c=getc(ficpar))=='#' && c!= EOF){ | while((c=getc(ficpar))=='#' && c!= EOF){ |
| Line 8054 Please run with mle=-1 to get a correct | Line 8478 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\ | 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); | Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
| }else | }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,\ | printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
| model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ | model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \ |
| jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); | jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2); |
| /*------------ free_vector -------------*/ | /*------------ free_vector -------------*/ |
| /* chdir(path); */ | /* chdir(path); */ |
| free_ivector(wav,1,imx); | /* free_ivector(wav,1,imx); */ /* Moved after last prevalence call */ |
| free_imatrix(dh,1,lastpass-firstpass+1,1,imx); | /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ |
| free_imatrix(bh,1,lastpass-firstpass+1,1,imx); | /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ |
| free_imatrix(mw,1,lastpass-firstpass+1,1,imx); | /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ |
| free_lvector(num,1,n); | free_lvector(num,1,n); |
| free_vector(agedc,1,n); | free_vector(agedc,1,n); |
| /*free_matrix(covar,0,NCOVMAX,1,n);*/ | /*free_matrix(covar,0,NCOVMAX,1,n);*/ |
| Line 8126 Please run with mle=-1 to get a correct | Line 8550 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",\ | /* 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); | 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) { | if (mobilav!=0) { |
| mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); | mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| Line 8385 Please run with mle=-1 to get a correct | Line 8814 Please run with mle=-1 to get a correct |
| if((nberr >0) || (nbwarn>0)){ | if((nberr >0) || (nbwarn>0)){ |
| printf("End of Imach with %d errors and/or %d warnings\n",nberr,nbwarn); | printf("End of Imach with %d errors and/or %d warnings. Please look at the log file for details.\n",nberr,nbwarn); |
| fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn); | fprintf(ficlog,"End of Imach with %d errors and/or warnings %d. Please look at the log file for details.\n",nberr,nbwarn); |
| }else{ | }else{ |
| printf("End of Imach\n"); | printf("End of Imach\n"); |
| fprintf(ficlog,"End of Imach\n"); | fprintf(ficlog,"End of Imach\n"); |