--- imach/src/imach.c 2016/08/26 14:23:35 1.238 +++ imach/src/imach.c 2016/09/07 14:10:18 1.248 @@ -1,6 +1,38 @@ -/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 brouard Exp $ +/* $Id: imach.c,v 1.248 2016/09/07 14:10:18 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.248 2016/09/07 14:10:18 brouard + *** empty log message *** + + Revision 1.247 2016/09/02 11:11:21 brouard + *** empty log message *** + + Revision 1.246 2016/09/02 08:49:22 brouard + *** empty log message *** + + Revision 1.245 2016/09/02 07:25:01 brouard + *** empty log message *** + + Revision 1.244 2016/09/02 07:17:34 brouard + *** empty log message *** + + Revision 1.243 2016/09/02 06:45:35 brouard + *** empty log message *** + + Revision 1.242 2016/08/30 15:01:20 brouard + Summary: Fixing a lots + + Revision 1.241 2016/08/29 17:17:25 brouard + Summary: gnuplot problem in Back projection to fix + + Revision 1.240 2016/08/29 07:53:18 brouard + Summary: Better + + Revision 1.239 2016/08/26 15:51:03 brouard + Summary: Improvement in Powell output in order to copy and paste + + Author: + Revision 1.238 2016/08/26 14:23:35 brouard Summary: Starting tests of 0.99 @@ -914,12 +946,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 brouard Exp $ */ +/* $Id: imach.c,v 1.248 2016/09/07 14:10:18 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; -char fullversion[]="$Revision: 1.238 $ $Date: 2016/08/26 14:23:35 $"; +char fullversion[]="$Revision: 1.248 $ $Date: 2016/09/07 14:10:18 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1032,7 +1064,8 @@ double dval; #define FTOL 1.0e-10 #define NRANSI -#define ITMAX 200 +#define ITMAX 200 +#define ITPOWMAX 20 /* This is now multiplied by the number of parameters */ #define TOL 2.0e-4 @@ -2061,10 +2094,10 @@ void powell(double p[], double **xi, int void linmin(double p[], double xi[], int n, double *fret, double (*func)(double [])); #else - void linmin(double p[], double xi[], int n, double *fret, - double (*func)(double []),int *flat); + void linmin(double p[], double xi[], int n, double *fret, + double (*func)(double []),int *flat); #endif - int i,ibig,j; + int i,ibig,j,jk,k; double del,t,*pt,*ptt,*xit; double directest; double fp,fptt; @@ -2096,31 +2129,64 @@ void powell(double p[], double **xi, int fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); /* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ for (i=1;i<=n;i++) { - printf(" %d %.12f",i, p[i]); - fprintf(ficlog," %d %.12lf",i, p[i]); fprintf(ficrespow," %.12lf", p[i]); } + fprintf(ficrespow,"\n");fflush(ficrespow); + printf("\n#model= 1 + age "); + fprintf(ficlog,"\n#model= 1 + age "); + if(nagesqr==1){ + printf(" + age*age "); + fprintf(ficlog," + age*age "); + } + for(j=1;j <=ncovmodel-2;j++){ + if(Typevar[j]==0) { + printf(" + V%d ",Tvar[j]); + fprintf(ficlog," + V%d ",Tvar[j]); + }else if(Typevar[j]==1) { + printf(" + V%d*age ",Tvar[j]); + fprintf(ficlog," + V%d*age ",Tvar[j]); + }else if(Typevar[j]==2) { + printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + } + } printf("\n"); +/* printf("12 47.0114589 0.0154322 33.2424412 0.3279905 2.3731903 */ +/* 13 -21.5392400 0.1118147 1.2680506 1.2973408 -1.0663662 */ fprintf(ficlog,"\n"); - fprintf(ficrespow,"\n");fflush(ficrespow); - if(*iter <=3){ + for(i=1,jk=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + for(j=1; j <=ncovmodel; j++){ + printf("%12.7f ",p[jk]); + fprintf(ficlog,"%12.7f ",p[jk]); + jk++; + } + printf("\n"); + fprintf(ficlog,"\n"); + } + } + } + if(*iter <=3 && *iter >1){ tml = *localtime(&rcurr_time); strcpy(strcurr,asctime(&tml)); rforecast_time=rcurr_time; itmp = strlen(strcurr); if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ - strcurr[itmp-1]='\0'; + strcurr[itmp-1]='\0'; printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); for(niterf=10;niterf<=30;niterf+=10){ - rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); - forecast_time = *localtime(&rforecast_time); - strcpy(strfor,asctime(&forecast_time)); - itmp = strlen(strfor); - if(strfor[itmp-1]=='\n') - strfor[itmp-1]='\0'; - printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); - fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); + rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); + forecast_time = *localtime(&rforecast_time); + strcpy(strfor,asctime(&forecast_time)); + itmp = strlen(strfor); + if(strfor[itmp-1]=='\n') + strfor[itmp-1]='\0'; + printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); + fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); } } for (i=1;i<=n;i++) { /* For each direction i */ @@ -2174,7 +2240,8 @@ void powell(double p[], double **xi, int /* printf("\n"); */ /* fprintf(ficlog,"\n"); */ } - if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */ + /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */ + if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */ /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */ /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */ @@ -2220,7 +2287,7 @@ void powell(double p[], double **xi, int free_vector(pt,1,n); return; } /* enough precision */ - if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); + if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations."); for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ ptt[j]=2.0*p[j]-pt[j]; xit[j]=p[j]-pt[j]; @@ -2500,7 +2567,7 @@ Earliest age to start was %d-%d=%d, ncvl /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */ /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */ - double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij) + double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij, int nres) { /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit matrix by transitions matrix until convergence is reached with precision ftolpl */ @@ -2523,6 +2590,7 @@ Earliest age to start was %d-%d=%d, ncvl /* If we start from prlim again, prlim tends to a constant matrix */ int i, ii,j,k; + int first=0; double *min, *max, *meandiff, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ double **out, cov[NCOVMAX+1], **bmij(); @@ -2559,15 +2627,49 @@ Earliest age to start was %d-%d=%d, ncvl cov[2]=agefin; if(nagesqr==1) cov[3]= agefin*agefin;; - for (k=1; k<=cptcovn;k++) { - /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; - /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */ + for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ + /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; + /* printf("bprevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ + } + /* for (k=1; k<=cptcovn;k++) { */ + /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ + /* /\* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); *\/ */ + /* } */ + for (k=1; k<=nsq;k++) { /* For single varying covariates only */ + /* Here comes the value of quantitative after renumbering k with single quantitative covariates */ + cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; + /* printf("prevalim Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ + } + /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; */ + /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */ + /* /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ + for (k=1; k<=cptcovage;k++){ /* For product with age */ + if(Dummy[Tvar[Tage[k]]]){ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + } else{ + cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; + } + /* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ + } + for (k=1; k<=cptcovprod;k++){ /* For product without age */ + /* printf("prevalim Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */ + if(Dummy[Tvard[k][1]==0]){ + if(Dummy[Tvard[k][2]==0]){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; + }else{ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; + } + }else{ + if(Dummy[Tvard[k][2]==0]){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; + }else{ + cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; + } + } } - for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; - for (k=1; k<=cptcovprod;k++) /* Useless */ - /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ @@ -2614,7 +2716,12 @@ Earliest age to start was %d-%d=%d, ncvl } } /* age loop */ /* After some age loop it doesn't converge */ - printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ + if(first){ + first=1; + printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\ +Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + } + fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ free_vector(min,1,nlstate); @@ -2747,7 +2854,8 @@ double **pmij(double **ps, double *cov, /* }else */ doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); }else{ - printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); + ; + /* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */ } } /*End ii */ } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ @@ -3130,7 +3238,7 @@ double func( double *x) Then computes with function pmij which return a matrix p[i][j] giving the elementary probability to be observed in j being in i according to the model. */ - ioffset=2+nagesqr+cptcovage; + ioffset=2+nagesqr ; /* Fixed */ for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */ cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ @@ -3149,7 +3257,8 @@ double func( double *x) */ for(mi=1; mi<= wav[i]-1; mi++){ for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ - cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; + /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */ + cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; } for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ @@ -3163,7 +3272,10 @@ double func( double *x) if(nagesqr==1) cov[3]= agexact*agexact; /* Should be changed here */ for (kk=1; kk<=cptcovage;kk++) { + if(!FixedV[Tvar[Tage[kk]]]) cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ + else + cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; } out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); @@ -3445,7 +3557,8 @@ double funcone( double *x) for(k=1; k<=nlstate; k++) ll[k]=0.; ioffset=0; for (i=1,ipmx=0, sw=0.; i<=imx; i++){ - ioffset=2+nagesqr+cptcovage; + /* ioffset=2+nagesqr+cptcovage; */ + ioffset=2+nagesqr; /* Fixed */ /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */ @@ -3472,46 +3585,51 @@ double funcone( double *x) for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ /* Wave varying (but not age varying) */ for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ - cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; - } + /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */ + cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; + } /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */ - /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ - /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ - /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */ - /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */ - /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ + /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ + /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ + /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */ + /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */ + /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */ - /* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */ - /* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */ - /* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */ + /* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */ + /* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */ + /* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */ /* } */ for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } 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 nlstate && (mle <5) ){ /* Jackson */ - 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 += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; - lli= log(survp); + for (j=1,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); }else if (mle==1){ - 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 */ } else if(mle==2){ - 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 */ } else if(mle==3){ /* exponential inter-extrapolation */ - lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ + lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ } else if (mle==4){ /* mle=4 no inter-extrapolation */ - lli=log(out[s1][s2]); /* Original formula */ + lli=log(out[s1][s2]); /* Original formula */ } else{ /* mle=0 back to 1 */ - lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ - /*lli=log(out[s1][s2]); */ /* Original formula */ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + /*lli=log(out[s1][s2]); */ /* Original formula */ } /* End of if */ ipmx +=1; sw += weight[i]; 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]); */ if(globpr){ - fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ + fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ - 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]); - for(k=1,llt=0.,l=0.; k<=nlstate; k++){ - llt +=ll[k]*gipmx/gsw; - fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); - } - fprintf(ficresilk," %10.6f\n", -llt); + 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]); + for(k=1,llt=0.,l=0.; k<=nlstate; k++){ + llt +=ll[k]*gipmx/gsw; + fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); + } + fprintf(ficresilk," %10.6f\n", -llt); } } /* end of wave */ } /* end of individual */ @@ -3905,8 +4023,8 @@ double hessij( double x[], double **hess kmax=kmax+10; } if(kmax >=10 || firstime ==1){ - printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol); - fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol); + printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); + fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); } @@ -4094,7 +4212,7 @@ void freqsummary(char fileres[], int ia 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)); @@ -4104,55 +4222,54 @@ void freqsummary(char fileres[], int ia } else{ fprintf(ficresphtm,"\nIMaCh PHTM_ %s\n %s
%s
\ -
\n\ +
\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 and dummy covariate value at beginning 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{ + } else{ fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s
\ -
\n\ +
\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); - + fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate)

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); j1=0; /* j=ncoveff; /\* Only fixed dummy covariates *\/ */ j=cptcoveff; /* Only dummy covariates of the model */ if (cptcovn<1) {j=1;ncodemax[1]=1;} - + first=1; - + /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: reference=low_education V1=0,V2=0 med_educ V1=1 V2=0, high_educ V1=0 V2=1 Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff */ - + for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */ posproptt=0.; /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ for (i=-5; i<=nlstate+ndeath; i++) for (jk=-5; jk<=nlstate+ndeath; jk++) - for(m=iagemin; m <= iagemax+3; m++) - freq[i][jk][m]=0; - + for(m=iagemin; m <= iagemax+3; m++) + freq[i][jk][m]=0; + for (i=1; i<=nlstate; i++) { for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0; + prop[i][m]=0; posprop[i]=0; pospropt[i]=0; } @@ -4162,7 +4279,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* meanqt[m][z1]=0.; */ /* } */ /* } */ - + dateintsum=0; k2cpt=0; /* For that combination of covariate j1, we count and print the frequencies in one pass */ @@ -4241,35 +4358,41 @@ Title=%s
Datafile=%s Firstpass=%d La } /* end iind = 1 to imx */ /* prop[s][age] is feeded for any initial and valid live state as well as freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ - - + + /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ pstamp(ficresp); - /* if (ncoveff>0) { */ - if (cptcoveff>0) { + if (cptcoveff>0){ fprintf(ficresp, "\n#********** Variable "); fprintf(ficresphtm, "\n

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

********** Variable "); + fprintf(ficlog, "\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)]); + if(DummyV[z1]){ + fprintf(ficresp, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficlog, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + }else{ + fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficlog, "V%d(varying)=%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(ficresphtm,""); 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, "",i,i); } fprintf(ficresp, "\n"); fprintf(ficresphtm, "\n"); - + /* Header of frequency table by age */ fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); fprintf(ficresphtmfr," "); @@ -4280,115 +4403,115 @@ Title=%s
Datafile=%s Firstpass=%d La } } fprintf(ficresphtmfr, "\n"); - + /* For each age */ for(iage=iagemin; iage <= iagemax+3; iage++){ fprintf(ficresphtm,""); if(iage==iagemax+1){ - fprintf(ficlog,"1"); - fprintf(ficresphtmfr," "); + fprintf(ficlog,"1"); + fprintf(ficresphtmfr," "); }else if(iage==iagemax+2){ - fprintf(ficlog,"0"); - fprintf(ficresphtmfr," "); + fprintf(ficlog,"0"); + fprintf(ficresphtmfr," "); }else if(iage==iagemax+3){ - fprintf(ficlog,"Total"); - fprintf(ficresphtmfr," "); + fprintf(ficlog,"Total"); + fprintf(ficresphtmfr," "); }else{ - if(first==1){ - first=0; - printf("See log file for details...\n"); - } - fprintf(ficresphtmfr," ",iage); - fprintf(ficlog,"Age %d", iage); + if(first==1){ + first=0; + printf("See log file for details...\n"); + } + fprintf(ficresphtmfr," ",iage); + fprintf(ficlog,"Age %d", iage); } for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][iage]; + for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) + pp[jk] += freq[jk][m][iage]; } for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][iage]; - if(pp[jk]>=1.e-10){ - if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); - } - fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); - }else{ - if(first==1) - printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - } + for(m=-1, pos=0; m <=0 ; m++) + pos += freq[jk][m][iage]; + if(pp[jk]>=1.e-10){ + if(first==1){ + printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + } + fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + }else{ + if(first==1) + printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); + fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); + } } - + for(jk=1; jk <=nlstate ; jk++){ - /* posprop[jk]=0; */ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ - pp[jk] += freq[jk][m][iage]; + /* posprop[jk]=0; */ + for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ + pp[jk] += freq[jk][m][iage]; } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ - + for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ - pos += pp[jk]; /* pos is the total number of transitions until this age */ - posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ - pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pos += pp[jk]; /* pos is the total number of transitions until this age */ + posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state + from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state + from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ } for(jk=1; jk <=nlstate ; jk++){ - if(pos>=1.e-5){ - if(first==1) - printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - }else{ - if(first==1) - printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - } - if( iage <= iagemax){ - if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - /*probs[iage][jk][j1]= pp[jk]/pos;*/ - /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ - } - else{ - fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage, prop[jk][iage],pospropta); - } - } - pospropt[jk] +=posprop[jk]; + if(pos>=1.e-5){ + if(first==1) + printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + }else{ + if(first==1) + printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); + fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); + } + if( iage <= iagemax){ + if(pos>=1.e-5){ + fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); + fprintf(ficresphtm,"",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); + /*probs[iage][jk][j1]= pp[jk]/pos;*/ + /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ + } + else{ + fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); + fprintf(ficresphtm,"",iage, prop[jk][iage],pospropta); + } + } + pospropt[jk] +=posprop[jk]; } /* end loop jk */ /* pospropt=0.; */ for(jk=-1; jk <=nlstate+ndeath; jk++){ - for(m=-1; m <=nlstate+ndeath; m++){ - if(freq[jk][m][iage] !=0 ) { /* minimizing output */ - if(first==1){ - printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); - } - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); - } - if(jk!=0 && m!=0) - fprintf(ficresphtmfr," ",freq[jk][m][iage]); - } + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[jk][m][iage] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); + } + fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); + } + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",freq[jk][m][iage]); + } } /* end loop jk */ posproptt=0.; for(jk=1; jk <=nlstate; jk++){ - posproptt += pospropt[jk]; + posproptt += pospropt[jk]; } fprintf(ficresphtmfr,"\n "); if(iage <= iagemax){ - fprintf(ficresp,"\n"); - fprintf(ficresphtm,"\n"); + fprintf(ficresp,"\n"); + fprintf(ficresphtm,"\n"); } if(first==1) - printf("Others in log...\n"); + printf("Others in log...\n"); fprintf(ficlog,"\n"); } /* end loop age iage */ fprintf(ficresphtm,""); for(jk=1; jk <=nlstate ; jk++){ if(posproptt < 1.e-5){ - fprintf(ficresphtm,"",pospropt[jk],posproptt); + fprintf(ficresphtm,"",pospropt[jk],posproptt); }else{ - fprintf(ficresphtm,"",pospropt[jk]/posproptt,pospropt[jk],posproptt); + fprintf(ficresphtm,"",pospropt[jk]/posproptt,pospropt[jk],posproptt); } } fprintf(ficresphtm,"\n"); @@ -4406,7 +4529,7 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtmfr,"
Age
0
0
Unknown
Unknown
Total
Total
%d
%d%d%.5f%.0f%.0f%dNaNq%.0f%.0f%d%.5f%.0f%.0f%dNaNq%.0f%.0f%.0f%.0f
TotNanq%.0f%.0fNanq%.0f%.0f%.5f%.0f%.0f%.5f%.0f%.0f
\n"); } /* end selected combination of covariate j1 */ dateintmean=dateintsum/k2cpt; - + fclose(ficresp); fclose(ficresphtm); fclose(ficresphtmfr); @@ -4741,173 +4864,173 @@ void concatwav(int wav[], int **dh, int /*********** Tricode ****************************/ void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum) -{ - /**< Uses cptcovn+2*cptcovprod as the number of covariates */ - /* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 - * Boring subroutine which should only output nbcode[Tvar[j]][k] - * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable - * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually); - */ + { + /**< Uses cptcovn+2*cptcovprod as the number of covariates */ + /* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 + * Boring subroutine which should only output nbcode[Tvar[j]][k] + * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable + * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually); + */ - int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX; - int modmaxcovj=0; /* Modality max of covariates j */ - int cptcode=0; /* Modality max of covariates j */ - int modmincovj=0; /* Modality min of covariates j */ + int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX; + int modmaxcovj=0; /* Modality max of covariates j */ + int cptcode=0; /* Modality max of covariates j */ + int modmincovj=0; /* Modality min of covariates j */ - /* cptcoveff=0; */ - /* *cptcov=0; */ + /* cptcoveff=0; */ + /* *cptcov=0; */ - for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ + for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ - /* Loop on covariates without age and products and no quantitative variable */ - /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */ - for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ - for (j=-1; (j < maxncov); j++) Ndum[j]=0; - if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ - switch(Fixed[k]) { - case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ - for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ - ij=(int)(covar[Tvar[k]][i]); - /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i - * If product of Vn*Vm, still boolean *: - * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables - * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */ - /* Finds for covariate j, n=Tvar[j] of Vn . ij is the - modality of the nth covariate of individual i. */ - if (ij > modmaxcovj) - modmaxcovj=ij; - else if (ij < modmincovj) - modmincovj=ij; - if ((ij < -1) && (ij > NCOVMAX)){ - printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); - exit(1); - }else - Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/ - /* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */ - /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ - /* getting the maximum value of the modality of the covariate - (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and - female ies 1, then modmaxcovj=1. - */ - } /* end for loop on individuals i */ - printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); - fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); - cptcode=modmaxcovj; - /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ - /*for (i=0; i<=cptcode; i++) {*/ - for (j=modmincovj; j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */ - printf("Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); - fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); - if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */ - if( j != -1){ - ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th - covariate for which somebody answered excluding - undefined. Usually 2: 0 and 1. */ - } - ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th - covariate for which somebody answered including - undefined. Usually 3: -1, 0 and 1. */ - } /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for - * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ - } /* Ndum[-1] number of undefined modalities */ + /* Loop on covariates without age and products and no quantitative variable */ + /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */ + for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ + for (j=-1; (j < maxncov); j++) Ndum[j]=0; + if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ + switch(Fixed[k]) { + case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ + for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ + ij=(int)(covar[Tvar[k]][i]); + /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i + * If product of Vn*Vm, still boolean *: + * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables + * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */ + /* Finds for covariate j, n=Tvar[j] of Vn . ij is the + modality of the nth covariate of individual i. */ + if (ij > modmaxcovj) + modmaxcovj=ij; + else if (ij < modmincovj) + modmincovj=ij; + if ((ij < -1) && (ij > NCOVMAX)){ + printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); + exit(1); + }else + Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/ + /* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */ + /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ + /* getting the maximum value of the modality of the covariate + (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and + female ies 1, then modmaxcovj=1. + */ + } /* end for loop on individuals i */ + printf(" Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); + fprintf(ficlog," Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); + cptcode=modmaxcovj; + /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ + /*for (i=0; i<=cptcode; i++) {*/ + for (j=modmincovj; j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */ + printf("Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); + fprintf(ficlog, "Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); + if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */ + if( j != -1){ + ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th + covariate for which somebody answered excluding + undefined. Usually 2: 0 and 1. */ + } + ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th + covariate for which somebody answered including + undefined. Usually 3: -1, 0 and 1. */ + } /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for + * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ + } /* Ndum[-1] number of undefined modalities */ - /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ - /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */ - /* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */ - /* modmincovj=3; modmaxcovj = 7; */ - /* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */ - /* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */ - /* defining two dummy variables: variables V1_1 and V1_2.*/ - /* nbcode[Tvar[j]][ij]=k; */ - /* nbcode[Tvar[j]][1]=0; */ - /* nbcode[Tvar[j]][2]=1; */ - /* nbcode[Tvar[j]][3]=2; */ - /* To be continued (not working yet). */ - ij=0; /* ij is similar to i but can jump over null modalities */ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ - if (Ndum[i] == 0) { /* If nobody responded to this modality k */ - break; - } - ij++; - nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/ - cptcode = ij; /* New max modality for covar j */ - } /* end of loop on modality i=-1 to 1 or more */ - break; - case 1: /* Testing on varying covariate, could be simple and - * should look at waves or product of fixed * - * varying. No time to test -1, assuming 0 and 1 only */ - ij=0; - for(i=0; i<=1;i++){ - nbcode[Tvar[k]][++ij]=i; - } - break; - default: - break; - } /* end switch */ - } /* end dummy test */ - - /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ - /* /\*recode from 0 *\/ */ - /* k is a modality. If we have model=V1+V1*sex */ - /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ - /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ - /* } */ - /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ - /* if (ij > ncodemax[j]) { */ - /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ - /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ - /* break; */ - /* } */ - /* } /\* end of loop on modality k *\/ */ - } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ - - for (k=-1; k< maxncov; k++) Ndum[k]=0; - /* Look at fixed dummy (single or product) covariates to check empty modalities */ - for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ - /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ - ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */ - Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */ - /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, {2, 1, 1, 1, 2, 1, 1, 0, 0} */ - } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */ - - ij=0; - /* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */ - for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ - /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ - /* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */ - if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */ - /* If product not in single variable we don't print results */ - /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ - ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */ - Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/ - Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */ - TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ - if(Fixed[k]!=0) - anyvaryingduminmodel=1; - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ - /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ - /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ - /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ - } - } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ - /* ij--; */ - /* cptcoveff=ij; /\*Number of total covariates*\/ */ - *cptcov=ij; /*Number of total real effective covariates: effective - * because they can be excluded from the model and real - * if in the model but excluded because missing values, but how to get k from ij?*/ - for(j=ij+1; j<= cptcovt; j++){ - Tvaraff[j]=0; - Tmodelind[j]=0; - } - for(j=ntveff+1; j<= cptcovt; j++){ - TmodelInvind[j]=0; - } - /* To be sorted */ - ; -} + /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ + /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */ + /* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */ + /* modmincovj=3; modmaxcovj = 7; */ + /* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */ + /* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */ + /* defining two dummy variables: variables V1_1 and V1_2.*/ + /* nbcode[Tvar[j]][ij]=k; */ + /* nbcode[Tvar[j]][1]=0; */ + /* nbcode[Tvar[j]][2]=1; */ + /* nbcode[Tvar[j]][3]=2; */ + /* To be continued (not working yet). */ + ij=0; /* ij is similar to i but can jump over null modalities */ + for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + if (Ndum[i] == 0) { /* If nobody responded to this modality k */ + break; + } + ij++; + nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/ + cptcode = ij; /* New max modality for covar j */ + } /* end of loop on modality i=-1 to 1 or more */ + break; + case 1: /* Testing on varying covariate, could be simple and + * should look at waves or product of fixed * + * varying. No time to test -1, assuming 0 and 1 only */ + ij=0; + for(i=0; i<=1;i++){ + nbcode[Tvar[k]][++ij]=i; + } + break; + default: + break; + } /* end switch */ + } /* end dummy test */ + + /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ + /* /\*recode from 0 *\/ */ + /* k is a modality. If we have model=V1+V1*sex */ + /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ + /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ + /* } */ + /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ + /* if (ij > ncodemax[j]) { */ + /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ + /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ + /* break; */ + /* } */ + /* } /\* end of loop on modality k *\/ */ + } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ + + for (k=-1; k< maxncov; k++) Ndum[k]=0; + /* Look at fixed dummy (single or product) covariates to check empty modalities */ + for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ + /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ + ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */ + Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */ + /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, {2, 1, 1, 1, 2, 1, 1, 0, 0} */ + } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */ + + ij=0; + /* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */ + for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ + /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ + /* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */ + if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */ + /* If product not in single variable we don't print results */ + /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ + ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */ + Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/ + Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */ + TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ + if(Fixed[k]!=0) + anyvaryingduminmodel=1; + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ + /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ + /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ + /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ + } + } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ + /* ij--; */ + /* cptcoveff=ij; /\*Number of total covariates*\/ */ + *cptcov=ij; /*Number of total real effective covariates: effective + * because they can be excluded from the model and real + * if in the model but excluded because missing values, but how to get k from ij?*/ + for(j=ij+1; j<= cptcovt; j++){ + Tvaraff[j]=0; + Tmodelind[j]=0; + } + for(j=ntveff+1; j<= cptcovt; j++){ + TmodelInvind[j]=0; + } + /* To be sorted */ + ; + } /*********** Health Expectancies ****************/ @@ -5364,7 +5487,7 @@ void concatwav(int wav[], int **dh, int xp[i] = x[i] + (i==theta ?delti[theta]:0); } - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nresult); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres); if (popbased==1) { if(mobilav ==0){ @@ -5396,7 +5519,7 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nresult); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres); if (popbased==1) { if(mobilav ==0){ @@ -5473,7 +5596,7 @@ void concatwav(int wav[], int **dh, int /* end ppptj */ /* x centered again */ - prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nresult); + prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nres); if (popbased==1) { if(mobilav ==0){ @@ -5568,7 +5691,9 @@ void concatwav(int wav[], int **dh, int pstamp(ficresvpl); fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); - fprintf(ficresvpl,"# Age"); + fprintf(ficresvpl,"# Age "); + if(nresult >=1) + fprintf(ficresvpl," Result# "); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); fprintf(ficresvpl,"\n"); @@ -5654,6 +5779,8 @@ void concatwav(int wav[], int **dh, int varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ fprintf(ficresvpl,"%.0f ",age ); + if(nresult >=1) + fprintf(ficresvpl,"%d ",nres ); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); fprintf(ficresvpl,"\n"); @@ -6058,7 +6185,7 @@ void printinghtml(char fileresu[], char jj1=0; for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k1=1; k1<=m;k1++){ + for(k1=1; k1<=m;k1++){ /* For each combination of covariate */ if(TKresult[nres]!= k1) continue; @@ -6086,51 +6213,51 @@ void printinghtml(char fileresu[], char } } /* aij, bij */ - 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); + 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-%d.svg
\ +",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); /* 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
\ -",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); + fprintf(fichtm,"
\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s_%d-2-%d.svg
\ +",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); /* Quasi-incidences */ 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
\ -",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); +divided by h: hPij/h : %s_%d-3-%d.svg
\ +",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); /* Survival functions (period) in state j */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. %s%d_%d.svg
\ -", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1); + fprintf(fichtm,"
\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\ +", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); } /* State specific survival functions (period) */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
\n- Survival functions from state %d in each live state and total.\ Or probability to survive in various states (1 to %d) being in state %d at different ages. \ - %s%d_%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1); + %s_%d%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); } /* Period (stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\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); + fprintf(fichtm,"
\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-%d.svg
\ +", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); } if(backcast==1){ /* Period (stable) back prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Convergence to period (stable) back 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,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1); + fprintf(fichtm,"
\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\ +", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); } } 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); + 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-%d.svg
\ +", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); } } 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); + 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-%d.svg
\ +",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres); } /* } /\* end i1 *\/ */ }/* End k1 */ @@ -6189,7 +6316,7 @@ See page 'Matrix of variance-covariance jj1=0; - for(nres=1; nres <= nresult; nres++) /* For each resultline */ + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ for(k1=1; k1<=m;k1++){ if(TKresult[nres]!= k1) continue; @@ -6213,17 +6340,18 @@ See page 'Matrix of variance-covariance } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"\n
- Observed (cross-sectional) and period (incidence based) \ -prevalence (with 95%% confidence interval) in state (%d): %s_%d-%d.svg\n
\ -",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); +prevalence (with 95%% confidence interval) in state (%d): %s_%d-%d-%d.svg\n
\ +",cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres); } fprintf(fichtm,"\n
- Total life expectancy by age and \ health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ true period expectancies (those weighted with period prevalences are also\ drawn in addition to the population based expectancies computed using\ - observed and cahotic prevalences: %s_%d.svg\n
\ -",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); + observed and cahotic prevalences: %s_%d-%d.svg\n
\ +",subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres); /* } /\* end i1 *\/ */ }/* End k1 */ + }/* End nres */ fprintf(fichtm,""); fflush(fichtm); } @@ -6289,7 +6417,7 @@ void printinggnuplot(char fileresu[], ch if(TKresult[nres]!= k1) continue; /* We are interested in selected combination by the resultline */ - printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); + /* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */ fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ @@ -6298,37 +6426,34 @@ void printinggnuplot(char fileresu[], ch /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */ /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */ - printf(" V%d=%d ",Tvaraff[k],vlv); + /* printf(" V%d=%d ",Tvaraff[k],vlv); */ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); } for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ - printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + /* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */ fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); } - printf("\n#\n"); + /* printf("\n#\n"); */ fprintf(ficgp,"\n#\n"); if(invalidvarcomb[k1]){ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); continue; } - 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\ -set ylabel \"Probability\" \n \ -set ter svg size 640, 480\n \ -plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1); + fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres); + fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); + fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); + fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); @@ -6336,9 +6461,9 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1)); if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ - fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */ + fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */ if(cptcoveff ==0){ - fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ", 2+(cpt-1), cpt ); + fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line lt 3", 2+(cpt-1), cpt ); }else{ kl=0; for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ @@ -6353,8 +6478,8 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u /*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){ - fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ - 4+(cpt-1), cpt ); /* 4 or 6 ?*/ + fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' w l lt 3",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ + 2+cptcoveff*2+(cpt-1), cpt ); /* 4 or 6 ?*/ }else{ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]); kl++; @@ -6393,7 +6518,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u continue; } - fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres); for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ if(vpopbased==0) fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); @@ -6423,7 +6548,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); } /* state */ } /* 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-%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */ } /* end nres */ } /* k1 end 2 eme*/ @@ -6431,7 +6556,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u /*3eme*/ for (k1=1; k1<= m ; k1 ++){ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k) + if(TKresult[nres]!= k1) continue; for (cpt=1; cpt<= nlstate ; cpt ++) { @@ -6455,7 +6580,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u /* 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); + fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres); fprintf(ficgp,"set ter svg size 640, 480\n\ plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt); /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); @@ -6501,7 +6626,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u continue; } - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); + fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); k=3; @@ -6547,7 +6672,7 @@ set ter svg size 640, 480\nunset log y\n continue; } - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); + fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); k=3; @@ -6602,7 +6727,7 @@ set ter svg size 640, 480\nunset log y\n continue; } - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); + fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1,nres); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); k=3; /* Offset */ @@ -6648,7 +6773,7 @@ set ter svg size 640, 480\nunset log y\n continue; } - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); + fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); k=3; /* Offset */ @@ -6700,7 +6825,7 @@ set ter svg size 640, 480\nunset log y\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,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ @@ -6816,7 +6941,7 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); } fprintf(ficgp,"\n#\n"); - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng); + fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng,nres); fprintf(ficgp,"\nset ter svg size 640, 480 "); if (ng==1){ fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ @@ -7791,14 +7916,36 @@ int readdata(char datafile[], int firsto /*-------- data file ----------*/ FILE *fic; char dummy[]=" "; - int i=0, j=0, n=0, iv=0; + int i=0, j=0, n=0, iv=0, v; int lstra; int linei, month, year,iout; char line[MAXLINE], linetmp[MAXLINE]; char stra[MAXLINE], strb[MAXLINE]; char *stratrunc; + DummyV=ivector(1,NCOVMAX); /* 1 to 3 */ + FixedV=ivector(1,NCOVMAX); /* 1 to 3 */ + for(v=1; v <=ncovcol;v++){ + DummyV[v]=0; + FixedV[v]=0; + } + for(v=ncovcol+1; v <=ncovcol+nqv;v++){ + DummyV[v]=1; + FixedV[v]=0; + } + for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){ + DummyV[v]=0; + FixedV[v]=1; + } + for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){ + DummyV[v]=1; + FixedV[v]=1; + } + for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){ + printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); + fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); + } if((fic=fopen(datafile,"r"))==NULL) { printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout); @@ -8419,27 +8566,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); - - for(v=1; v <=ncovcol;v++){ - DummyV[v]=0; - FixedV[v]=0; - } - for(v=ncovcol+1; v <=ncovcol+nqv;v++){ - DummyV[v]=1; - FixedV[v]=0; - } - for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){ - DummyV[v]=0; - FixedV[v]=1; - } - for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){ - DummyV[v]=1; - FixedV[v]=1; - } - for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){ - printf("Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); - fprintf(ficlog,"Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); - } + for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ Fixed[k]= 0; @@ -8464,7 +8591,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( TvarFind[ncovf]=k; TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ - }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */ + }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */ Fixed[k]= 0; Dummy[k]= 1; nqfveff++; @@ -8478,7 +8605,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( TvarFind[ncovf]=k; TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ - }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying variables */ + }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */ Fixed[k]= 1; Dummy[k]= 0; ntveff++; /* Only simple time varying dummy variable */ @@ -8489,7 +8616,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( TvarsDind[nsd]=k; ncovv++; /* Only simple time varying variables */ TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; + TvarVind[ncovv]=k; /* TvarVind[2]=2 TvarVind[3]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv); @@ -8505,7 +8632,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( TvarsQ[nsq]=Tvar[k]; TvarsQind[nsq]=k; TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; + TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */ @@ -8517,167 +8644,167 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( TvarA[ncova]=Tvar[k]; TvarAind[ncova]=k; if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */ - Fixed[k]= 2; - Dummy[k]= 2; - modell[k].maintype= ATYPE; - modell[k].subtype= APFD; - /* ncoveff++; */ + Fixed[k]= 2; + Dummy[k]= 2; + modell[k].maintype= ATYPE; + modell[k].subtype= APFD; + /* ncoveff++; */ }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/ - Fixed[k]= 2; - Dummy[k]= 3; - modell[k].maintype= ATYPE; - modell[k].subtype= APFQ; /* Product age * fixed quantitative */ - /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */ + Fixed[k]= 2; + Dummy[k]= 3; + modell[k].maintype= ATYPE; + modell[k].subtype= APFQ; /* Product age * fixed quantitative */ + /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */ }else if( Tvar[k] <=ncovcol+nqv+ntv ){ - Fixed[k]= 3; - Dummy[k]= 2; - modell[k].maintype= ATYPE; - modell[k].subtype= APVD; /* Product age * varying dummy */ - /* ntveff++; /\* Only simple time varying dummy variable *\/ */ + Fixed[k]= 3; + Dummy[k]= 2; + modell[k].maintype= ATYPE; + modell[k].subtype= APVD; /* Product age * varying dummy */ + /* ntveff++; /\* Only simple time varying dummy variable *\/ */ }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 3; - Dummy[k]= 3; - modell[k].maintype= ATYPE; - modell[k].subtype= APVQ; /* Product age * varying quantitative */ - /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */ + Fixed[k]= 3; + Dummy[k]= 3; + modell[k].maintype= ATYPE; + modell[k].subtype= APVQ; /* Product age * varying quantitative */ + /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */ } }else if (Typevar[k] == 2) { /* product without age */ k1=Tposprod[k]; if(Tvard[k1][1] <=ncovcol){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 1; - Dummy[k]= 0; - modell[k].maintype= FTYPE; - modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */ - ncovf++; /* Fixed variables without age */ - TvarF[ncovf]=Tvar[k]; - TvarFind[ncovf]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[k]= 0; /* or 2 ?*/ - Dummy[k]= 1; - modell[k].maintype= FTYPE; - modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */ - ncovf++; /* Varying variables without age */ - TvarF[ncovf]=Tvar[k]; - TvarFind[ncovf]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 0; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */ + ncovf++; /* Fixed variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 0; /* or 2 ?*/ + Dummy[k]= 1; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */ + ncovf++; /* Varying variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } }else if(Tvard[k1][1] <=ncovcol+nqv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 0; /* or 2 ?*/ - Dummy[k]= 1; - modell[k].maintype= FTYPE; - modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */ - ncovf++; /* Fixed variables without age */ - TvarF[ncovf]=Tvar[k]; - TvarFind[ncovf]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 0; /* or 2 ?*/ + Dummy[k]= 1; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */ + ncovf++; /* Fixed variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 0; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - modell[k].maintype= VTYPE; - modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */ - ncovv++; /* Varying variables without age */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } }else{ - printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); - fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); - } /* end k1 */ + printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); + fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); + } /*end k1*/ }else{ printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); @@ -8719,7 +8846,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( } int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) -{ +{/* Check ages at death */ int i, m; int firstone=0; @@ -9192,14 +9319,14 @@ int back_prevalence_limit(double *p, dou if(mobilavproj > 0){ /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */ /* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */ - bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k); + bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k, nres); }else if (mobilavproj == 0){ printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj); fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj); exit(1); }else{ /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */ - bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k); + bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k,nres); } fprintf(ficresplb,"%.0f ",age ); for(j=1;j<=cptcoveff;j++) @@ -9763,6 +9890,12 @@ int main(int argc, char *argv[]) delti=delti3[1][1]; /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ if(mle==-1){ /* Print a wizard for help writing covariance matrix */ +/* We could also provide initial parameters values giving by simple logistic regression + * only one way, that is without matrix product. We will have nlstate maximizations */ + /* for(i=1;i=1){ /* Could be 1 or 2, Real Maximization */ /* mlikeli uses func not funcone */ + /* for(i=1;i