--- imach/src/imach.c 2016/02/15 23:22:02 1.220 +++ imach/src/imach.c 2016/02/17 08:14:50 1.222 @@ -1,8 +1,11 @@ -/* $Id: imach.c,v 1.220 2016/02/15 23:22:02 brouard Exp $ +/* $Id: imach.c,v 1.222 2016/02/17 08:14:50 brouard Exp $ $State: Exp $ $Log: imach.c,v $ - Revision 1.220 2016/02/15 23:22:02 brouard - Summary: 0.99r2 + Revision 1.222 2016/02/17 08:14:50 brouard + Summary: Probably last 0.98 stable version 0.98r6 + + Revision 1.221 2016/02/15 23:35:36 brouard + Summary: minor bug Revision 1.219 2016/02/15 00:48:12 brouard *** empty log message *** @@ -831,12 +834,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.220 2016/02/15 23:22:02 brouard Exp $ */ +/* $Id: imach.c,v 1.222 2016/02/17 08:14:50 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.220 $ $Date: 2016/02/15 23:22:02 $"; +char fullversion[]="$Revision: 1.222 $ $Date: 2016/02/17 08:14:50 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -2398,60 +2401,60 @@ double **pmij(double **ps, double *cov, /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) { - /* Computes the backward probability at age agefin and covariate ij - * and returns in **ps as well as **bmij. - */ + /* Computes the backward probability at age agefin and covariate ij + * and returns in **ps as well as **bmij. + */ int i, ii, j,k; - - double **out, **pmij(); - double sumnew=0.; + + double **out, **pmij(); + double sumnew=0.; double agefin; - - double **dnewm, **dsavm, **doldm; - double **bbmij; - + + double **dnewm, **dsavm, **doldm; + double **bbmij; + doldm=ddoldms; /* global pointers */ - dnewm=ddnewms; - dsavm=ddsavms; - - agefin=cov[2]; - /* bmij *//* age is cov[2], ij is included in cov, but we need for - the observed prevalence (with this covariate ij) */ - dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); - /* We do have the matrix Px in savm and we need pij */ - for (j=1;j<=nlstate+ndeath;j++){ - sumnew=0.; /* w1 p11 + w2 p21 only on live states */ - for (ii=1;ii<=nlstate;ii++){ - sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; - } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ - for (ii=1;ii<=nlstate+ndeath;ii++){ - if(sumnew >= 1.e-10){ - /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ - /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ - /* }else if(agefin >= agemaxpar+stepm/YEARM){ */ - /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ - /* }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); - } - } /*End ii */ - } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ - /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */ - bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ - /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ - /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ - /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ - /* left Product of this matrix by diag matrix of prevalences (savm) */ - for (j=1;j<=nlstate+ndeath;j++){ - for (ii=1;ii<=nlstate+ndeath;ii++){ - dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); - } - } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */ - ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ - /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ - /* end bmij */ - return ps; + dnewm=ddnewms; + dsavm=ddsavms; + + agefin=cov[2]; + /* bmij *//* age is cov[2], ij is included in cov, but we need for + the observed prevalence (with this covariate ij) */ + dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); + /* We do have the matrix Px in savm and we need pij */ + for (j=1;j<=nlstate+ndeath;j++){ + sumnew=0.; /* w1 p11 + w2 p21 only on live states */ + for (ii=1;ii<=nlstate;ii++){ + sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; + } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ + for (ii=1;ii<=nlstate+ndeath;ii++){ + if(sumnew >= 1.e-10){ + /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ + /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ + /* }else if(agefin >= agemaxpar+stepm/YEARM){ */ + /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ + /* }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); + } + } /*End ii */ + } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ + /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */ + bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ + /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ + /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ + /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ + /* left Product of this matrix by diag matrix of prevalences (savm) */ + for (j=1;j<=nlstate+ndeath;j++){ + for (ii=1;ii<=nlstate+ndeath;ii++){ + dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); + } + } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */ + ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ + /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ + /* end bmij */ + return ps; } /*************** transition probabilities ***************/ @@ -2657,7 +2660,7 @@ double ***hpxij(double ***po, int nhstep /************* Higher Back Matrix Product ***************/ /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */ - double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) +double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) { /* Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' months (i.e. until @@ -2669,16 +2672,16 @@ double ***hpxij(double ***po, int nhstep Model is determined by parameters x and covariates have to be included manually here. - */ + */ int i, j, d, h, k; double **out, cov[NCOVMAX+1]; double **newm; double agexact; double agebegin, ageend; - double **oldm, **savm; + double **oldm, **savm; - oldm=oldms;savm=savms; + oldm=oldms;savm=savms; /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) for (j=1;j<=nlstate+ndeath;j++){ @@ -2695,27 +2698,27 @@ double ***hpxij(double ***po, int nhstep /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ cov[2]=agexact; if(nagesqr==1) - cov[3]= agexact*agexact; + cov[3]= agexact*agexact; for (k=1; k<=cptcovn;k++) - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; - /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ - /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; - /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ + /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; - /* 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)]; + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ /* Careful transposed matrix */ - /* age is in cov[2] */ + /* age is in cov[2] */ /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ - /* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ + /* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ - 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); + 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* if((int)age == 70){ */ /* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */ /* for(i=1; i<=nlstate+ndeath; i++) { */ @@ -2735,12 +2738,12 @@ double ***hpxij(double ***po, int nhstep } for(i=1; i<=nlstate+ndeath; i++) for(j=1;j<=nlstate+ndeath;j++) { - po[i][j][h]=newm[i][j]; - /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/ + po[i][j][h]=newm[i][j]; + /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/ } /*printf("h=%d ",h);*/ } /* end h */ - /* printf("\n H=%d \n",h); */ + /* printf("\n H=%d \n",h); */ return po; } @@ -3962,9 +3965,9 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtm,"Tot"); for(jk=1; jk <=nlstate ; jk++){ if(posproptt < 1.e-5){ - fprintf(ficresphtm,"%.5f%.0f%.0f",pospropt[jk]/posproptt,pospropt[jk],posproptt); + fprintf(ficresphtm,"Nanq%.0f%.0f",pospropt[jk],posproptt); }else{ - fprintf(ficresphtm,"Nanq%.0f%.0f",pospropt[jk]/posproptt,pospropt[jk],posproptt); + fprintf(ficresphtm,"%.5f%.0f%.0f",pospropt[jk]/posproptt,pospropt[jk],posproptt); } } fprintf(ficresphtm,"\n"); @@ -3995,97 +3998,97 @@ Title=%s
Datafile=%s Firstpass=%d La } /************ Prevalence ********************/ -void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) -{ - /* 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. - */ + void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) + { + /* 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. + */ - int i, m, jk, j1, bool, z1,j; - int mi; /* Effective wave */ - int iage; - double agebegin, ageend; - - double **prop; - double posprop; - double y2; /* in fractional years */ - int iagemin, iagemax; - int first; /** to stop verbosity which is redirected to log file */ - - iagemin= (int) agemin; - iagemax= (int) agemax; - /*pp=vector(1,nlstate);*/ - prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); - /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ - j1=0; - - /*j=cptcoveff;*/ - if (cptcovn<1) {j=1;ncodemax[1]=1;} - - first=1; - for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ - for (i=1; i<=nlstate; i++) - for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) - prop[i][iage]=0.0; - - for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; - if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/ - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) - bool=0; - } - if (bool==1) { /* For this combination of covariates values, this individual fits */ - /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ - for(mi=1; mi=firstpass && m <=lastpass){ - y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ - if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ - if(agev[m][i]==0) agev[m][i]=iagemax+1; - if(agev[m][i]==1) agev[m][i]=iagemax+2; - if((int)agev[m][i] iagemax+3+AGEMARGE){ - printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); - exit(1); - } - if (s[m][i]>0 && s[m][i]<=nlstate) { - /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ - prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ - prop[s[m][i]][iagemax+3] += weight[i]; - } /* end valid statuses */ - } /* end selection of dates */ - } /* end selection of waves */ - } /* end effective waves */ - } /* end bool */ - } - for(i=iagemin; i <= iagemax+3; i++){ - for(jk=1,posprop=0; jk <=nlstate ; jk++) { - posprop += prop[jk][i]; - } + int i, m, jk, j1, bool, z1,j; + int mi; /* Effective wave */ + int iage; + double agebegin, ageend; + + double **prop; + double posprop; + double y2; /* in fractional years */ + int iagemin, iagemax; + int first; /** to stop verbosity which is redirected to log file */ + + iagemin= (int) agemin; + iagemax= (int) agemax; + /*pp=vector(1,nlstate);*/ + prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); + /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ + j1=0; + + /*j=cptcoveff;*/ + if (cptcovn<1) {j=1;ncodemax[1]=1;} + + first=1; + for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ + for (i=1; i<=nlstate; i++) + for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) + prop[i][iage]=0.0; + + for (i=1; i<=imx; i++) { /* Each individual */ + bool=1; + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/ + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) + bool=0; + } + if (bool==1) { /* For this combination of covariates values, this individual fits */ + /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ + for(mi=1; mi=firstpass && m <=lastpass){ + y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ + if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ + if(agev[m][i]==0) agev[m][i]=iagemax+1; + if(agev[m][i]==1) agev[m][i]=iagemax+2; + if((int)agev[m][i] iagemax+3+AGEMARGE){ + printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); + exit(1); + } + if (s[m][i]>0 && s[m][i]<=nlstate) { + /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ + prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ + prop[s[m][i]][iagemax+3] += weight[i]; + } /* end valid statuses */ + } /* end selection of dates */ + } /* end selection of waves */ + } /* end effective waves */ + } /* end bool */ + } + for(i=iagemin; i <= iagemax+3; i++){ + for(jk=1,posprop=0; jk <=nlstate ; jk++) { + posprop += prop[jk][i]; + } - for(jk=1; jk <=nlstate ; jk++){ - if( i <= iagemax){ - if(posprop>=1.e-5){ - probs[i][jk][j1]= prop[jk][i]/posprop; - } else{ - if(first==1){ - first=0; - printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); - } - } - } - }/* end jk */ - }/* end i */ - /*} *//* end i1 */ - } /* end j1 */ - - /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ - /*free_vector(pp,1,nlstate);*/ - free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE); -} /* End of prevalence */ + for(jk=1; jk <=nlstate ; jk++){ + if( i <= iagemax){ + if(posprop>=1.e-5){ + probs[i][jk][j1]= prop[jk][i]/posprop; + } else{ + if(first==1){ + first=0; + printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); + } + } + } + }/* end jk */ + }/* end i */ + /*} *//* end i1 */ + } /* end j1 */ + + /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ + /*free_vector(pp,1,nlstate);*/ + free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE); + } /* End of prevalence */ /************* Waves Concatenation ***************/ @@ -4520,7 +4523,7 @@ void cvevsij(double ***eij, double x[], { /* Covariances of health expectancies eij and of total life expectancies according - to initial status i, ei. . + to initial status i, ei. . */ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; int nhstepma, nstepma; /* Decreasing with age */ @@ -4626,47 +4629,47 @@ void cvevsij(double ***eij, double x[], decrease memory allocation */ for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++){ - xp[i] = x[i] + (i==theta ?delti[theta]:0); - xm[i] = x[i] - (i==theta ?delti[theta]:0); + xp[i] = x[i] + (i==theta ?delti[theta]:0); + xm[i] = x[i] - (i==theta ?delti[theta]:0); } hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij); hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij); for(j=1; j<= nlstate; j++){ - for(i=1; i<=nlstate; i++){ - for(h=0; h<=nhstepm-1; h++){ - gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.; - gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.; - } - } + for(i=1; i<=nlstate; i++){ + for(h=0; h<=nhstepm-1; h++){ + gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.; + gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.; + } + } } for(ij=1; ij<= nlstate*nlstate; ij++) - for(h=0; h<=nhstepm-1; h++){ - gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta]; - } + for(h=0; h<=nhstepm-1; h++){ + gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta]; + } }/* End theta */ for(h=0; h<=nhstepm-1; h++) for(j=1; j<=nlstate*nlstate;j++) - for(theta=1; theta <=npar; theta++) - trgradg[h][j][theta]=gradg[h][theta][j]; + for(theta=1; theta <=npar; theta++) + trgradg[h][j][theta]=gradg[h][theta][j]; - for(ij=1;ij<=nlstate*nlstate;ij++) + for(ij=1;ij<=nlstate*nlstate;ij++) for(ji=1;ji<=nlstate*nlstate;ji++) - varhe[ij][ji][(int)age] =0.; + varhe[ij][ji][(int)age] =0.; - printf("%d|",(int)age);fflush(stdout); - fprintf(ficlog,"%d|",(int)age);fflush(ficlog); - for(h=0;h<=nhstepm-1;h++){ + printf("%d|",(int)age);fflush(stdout); + fprintf(ficlog,"%d|",(int)age);fflush(ficlog); + for(h=0;h<=nhstepm-1;h++){ for(k=0;k<=nhstepm-1;k++){ - matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]); - for(ij=1;ij<=nlstate*nlstate;ij++) - for(ji=1;ji<=nlstate*nlstate;ji++) - varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; + matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]); + for(ij=1;ij<=nlstate*nlstate;ij++) + for(ji=1;ji<=nlstate*nlstate;ji++) + varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; } } @@ -4674,22 +4677,22 @@ void cvevsij(double ***eij, double x[], hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) - for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ - eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf; + for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ + eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf; - /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ - } + } fprintf(ficresstdeij,"%3.0f",age ); for(i=1; i<=nlstate;i++){ eip=0.; vip=0.; for(j=1; j<=nlstate;j++){ - eip += eij[i][j][(int)age]; - for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */ - vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age]; - fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) ); + eip += eij[i][j][(int)age]; + for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */ + vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age]; + fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) ); } fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip)); } @@ -4698,13 +4701,13 @@ void cvevsij(double ***eij, double x[], fprintf(ficrescveij,"%3.0f",age ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ - cptj= (j-1)*nlstate+i; - for(i2=1; i2<=nlstate;i2++) - for(j2=1; j2<=nlstate;j2++){ - cptj2= (j2-1)*nlstate+i2; - if(cptj2 <= cptj) - fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]); - } + cptj= (j-1)*nlstate+i; + for(i2=1; i2<=nlstate;i2++) + for(j2=1; j2<=nlstate;j2++){ + cptj2= (j2-1)*nlstate+i2; + if(cptj2 <= cptj) + fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]); + } } fprintf(ficrescveij,"\n"); @@ -5161,86 +5164,86 @@ void cvevsij(double ***eij, double x[], /************ Variance of one-step probabilities ******************/ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) -{ - int i, j=0, k1, l1, tj; - int k2, l2, j1, z1; - int k=0, l; - int first=1, first1, first2; - double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; - double **dnewm,**doldm; - double *xp; - double *gp, *gm; - double **gradg, **trgradg; - double **mu; - double age, cov[NCOVMAX+1]; - double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ - int theta; - char fileresprob[FILENAMELENGTH]; - char fileresprobcov[FILENAMELENGTH]; - char fileresprobcor[FILENAMELENGTH]; - double ***varpij; - - strcpy(fileresprob,"PROB_"); - strcat(fileresprob,fileres); - if((ficresprob=fopen(fileresprob,"w"))==NULL) { - printf("Problem with resultfile: %s\n", fileresprob); - fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); - } - strcpy(fileresprobcov,"PROBCOV_"); - strcat(fileresprobcov,fileresu); - if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { - printf("Problem with resultfile: %s\n", fileresprobcov); - fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); - } - strcpy(fileresprobcor,"PROBCOR_"); - strcat(fileresprobcor,fileresu); - if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { - printf("Problem with resultfile: %s\n", fileresprobcor); - fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); - } - printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); - fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); - printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); - fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); - printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); - fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); - pstamp(ficresprob); - fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); - fprintf(ficresprob,"# Age"); - pstamp(ficresprobcov); - fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); - fprintf(ficresprobcov,"# Age"); - pstamp(ficresprobcor); - fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); - fprintf(ficresprobcor,"# Age"); - + { + int i, j=0, k1, l1, tj; + int k2, l2, j1, z1; + int k=0, l; + int first=1, first1, first2; + double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; + double **dnewm,**doldm; + double *xp; + double *gp, *gm; + double **gradg, **trgradg; + double **mu; + double age, cov[NCOVMAX+1]; + double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ + int theta; + char fileresprob[FILENAMELENGTH]; + char fileresprobcov[FILENAMELENGTH]; + char fileresprobcor[FILENAMELENGTH]; + double ***varpij; + + strcpy(fileresprob,"PROB_"); + strcat(fileresprob,fileres); + if((ficresprob=fopen(fileresprob,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprob); + fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); + } + strcpy(fileresprobcov,"PROBCOV_"); + strcat(fileresprobcov,fileresu); + if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprobcov); + fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); + } + strcpy(fileresprobcor,"PROBCOR_"); + strcat(fileresprobcor,fileresu); + if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprobcor); + fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); + } + printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); + fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); + printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); + fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); + printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); + fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); + pstamp(ficresprob); + fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); + fprintf(ficresprob,"# Age"); + pstamp(ficresprobcov); + fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); + fprintf(ficresprobcov,"# Age"); + pstamp(ficresprobcor); + fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); + fprintf(ficresprobcor,"# Age"); - for(i=1; i<=nlstate;i++) - for(j=1; j<=(nlstate+ndeath);j++){ - fprintf(ficresprob," p%1d-%1d (SE)",i,j); - fprintf(ficresprobcov," p%1d-%1d ",i,j); - fprintf(ficresprobcor," p%1d-%1d ",i,j); - } - /* fprintf(ficresprob,"\n"); - fprintf(ficresprobcov,"\n"); - fprintf(ficresprobcor,"\n"); - */ - xp=vector(1,npar); - dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); - doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); - mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); - varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage); - first=1; - fprintf(ficgp,"\n# Routine varprob"); - fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); - fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.
    • \n",optionfilehtmcov); - fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); - fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ + for(i=1; i<=nlstate;i++) + for(j=1; j<=(nlstate+ndeath);j++){ + fprintf(ficresprob," p%1d-%1d (SE)",i,j); + fprintf(ficresprobcov," p%1d-%1d ",i,j); + fprintf(ficresprobcor," p%1d-%1d ",i,j); + } + /* fprintf(ficresprob,"\n"); + fprintf(ficresprobcov,"\n"); + fprintf(ficresprobcor,"\n"); + */ + xp=vector(1,npar); + dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); + doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); + mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); + varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage); + first=1; + fprintf(ficgp,"\n# Routine varprob"); + fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); + fprintf(fichtm,"\n"); + + fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.
    • \n",optionfilehtmcov); + fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); + fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ and drawn. It helps understanding how is the covariance between two incidences.\ They are expressed in year-1 in order to be less dependent of stepm.
    \n"); - fprintf(fichtmcov,"\n
    Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \ + fprintf(fichtmcov,"\n
    Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \ It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \ would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 \ standard deviations wide on each axis.
    \ @@ -5248,252 +5251,252 @@ standard deviations wide on each axis. < and made the appropriate rotation to look at the uncorrelated principal directions.
    \ To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.
    \n"); - cov[1]=1; - /* tj=cptcoveff; */ - tj = (int) pow(2,cptcoveff); - if (cptcovn<1) {tj=1;ncodemax[1]=1;} - j1=0; - for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates */ - if (cptcovn>0) { - fprintf(ficresprob, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresprob, "**********\n#\n"); - fprintf(ficresprobcov, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresprobcov, "**********\n#\n"); + cov[1]=1; + /* tj=cptcoveff; */ + tj = (int) pow(2,cptcoveff); + if (cptcovn<1) {tj=1;ncodemax[1]=1;} + j1=0; + for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates */ + if (cptcovn>0) { + fprintf(ficresprob, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresprob, "**********\n#\n"); + fprintf(ficresprobcov, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresprobcov, "**********\n#\n"); - fprintf(ficgp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficgp, "**********\n#\n"); + fprintf(ficgp, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficgp, "**********\n#\n"); - fprintf(fichtmcov, "\n
    ********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(fichtmcov, "**********\n
    "); + fprintf(fichtmcov, "\n
    ********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(fichtmcov, "**********\n
    "); - fprintf(ficresprobcor, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresprobcor, "**********\n#"); - if(invalidvarcomb[j1]){ - fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); - fprintf(fichtmcov,"\n

    Combination (%d) ignored because no cases

    \n",j1); - continue; - } - } - gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath)); - trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); - gp=vector(1,(nlstate)*(nlstate+ndeath)); - gm=vector(1,(nlstate)*(nlstate+ndeath)); - for (age=bage; age<=fage; age ++){ - cov[2]=age; - if(nagesqr==1) - cov[3]= age*age; - for (k=1; k<=cptcovn;k++) { - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; - /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 - * 1 1 1 1 1 - * 2 2 1 1 1 - * 3 1 2 1 1 - */ - /* nbcode[1][1]=0 nbcode[1][2]=1;*/ - } - /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; - for (k=1; k<=cptcovprod;k++) - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; + fprintf(ficresprobcor, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresprobcor, "**********\n#"); + if(invalidvarcomb[j1]){ + fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); + fprintf(fichtmcov,"\n

    Combination (%d) ignored because no cases

    \n",j1); + continue; + } + } + gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath)); + trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); + gp=vector(1,(nlstate)*(nlstate+ndeath)); + gm=vector(1,(nlstate)*(nlstate+ndeath)); + for (age=bage; age<=fage; age ++){ + cov[2]=age; + if(nagesqr==1) + cov[3]= age*age; + for (k=1; k<=cptcovn;k++) { + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; + /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 + * 1 1 1 1 1 + * 2 2 1 1 1 + * 3 1 2 1 1 + */ + /* nbcode[1][1]=0 nbcode[1][2]=1;*/ + } + /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + for (k=1; k<=cptcovprod;k++) + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; - for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++) - xp[i] = x[i] + (i==theta ?delti[theta]:(double)0); + for(theta=1; theta <=npar; theta++){ + for(i=1; i<=npar; i++) + xp[i] = x[i] + (i==theta ?delti[theta]:(double)0); - pmij(pmmij,cov,ncovmodel,xp,nlstate); + pmij(pmmij,cov,ncovmodel,xp,nlstate); - k=0; - for(i=1; i<= (nlstate); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gp[k]=pmmij[i][j]; - } - } + k=0; + for(i=1; i<= (nlstate); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gp[k]=pmmij[i][j]; + } + } - for(i=1; i<=npar; i++) - xp[i] = x[i] - (i==theta ?delti[theta]:(double)0); + for(i=1; i<=npar; i++) + xp[i] = x[i] - (i==theta ?delti[theta]:(double)0); - pmij(pmmij,cov,ncovmodel,xp,nlstate); - k=0; - for(i=1; i<=(nlstate); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gm[k]=pmmij[i][j]; - } - } + pmij(pmmij,cov,ncovmodel,xp,nlstate); + k=0; + for(i=1; i<=(nlstate); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gm[k]=pmmij[i][j]; + } + } - for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) - gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta]; - } + for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) + gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta]; + } - for(j=1; j<=(nlstate)*(nlstate+ndeath);j++) - for(theta=1; theta <=npar; theta++) - trgradg[j][theta]=gradg[theta][j]; + for(j=1; j<=(nlstate)*(nlstate+ndeath);j++) + for(theta=1; theta <=npar; theta++) + trgradg[j][theta]=gradg[theta][j]; - matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg); + matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg); - pmij(pmmij,cov,ncovmodel,x,nlstate); + pmij(pmmij,cov,ncovmodel,x,nlstate); - k=0; - for(i=1; i<=(nlstate); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - mu[k][(int) age]=pmmij[i][j]; - } - } - for(i=1;i<=(nlstate)*(nlstate+ndeath);i++) - for(j=1;j<=(nlstate)*(nlstate+ndeath);j++) - varpij[i][j][(int)age] = doldm[i][j]; + k=0; + for(i=1; i<=(nlstate); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + mu[k][(int) age]=pmmij[i][j]; + } + } + for(i=1;i<=(nlstate)*(nlstate+ndeath);i++) + for(j=1;j<=(nlstate)*(nlstate+ndeath);j++) + varpij[i][j][(int)age] = doldm[i][j]; - /*printf("\n%d ",(int)age); - for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ - printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); - fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); - }*/ + /*printf("\n%d ",(int)age); + for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ + printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); + fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); + }*/ - fprintf(ficresprob,"\n%d ",(int)age); - fprintf(ficresprobcov,"\n%d ",(int)age); - fprintf(ficresprobcor,"\n%d ",(int)age); + fprintf(ficresprob,"\n%d ",(int)age); + fprintf(ficresprobcov,"\n%d ",(int)age); + fprintf(ficresprobcor,"\n%d ",(int)age); - for (i=1; i<=(nlstate)*(nlstate+ndeath);i++) - fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age])); - for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ - fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]); - fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]); - } - i=0; - for (k=1; k<=(nlstate);k++){ - for (l=1; l<=(nlstate+ndeath);l++){ - i++; - fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l); - fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l); - for (j=1; j<=i;j++){ - /* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */ - fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]); - fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age])); - } - } - }/* end of loop for state */ - } /* end of loop for age */ - free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); - free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); - free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); - free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); - - /* Confidence intervalle of pij */ - /* - fprintf(ficgp,"\nunset parametric;unset label"); - fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\""); - fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); - fprintf(fichtm,"\n
    Probability with confidence intervals expressed in year-1 :pijgr%s.png, ",optionfilefiname,optionfilefiname); - fprintf(fichtm,"\n
    ",optionfilefiname); - fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname); - fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob); - */ - - /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ - first1=1;first2=2; - for (k2=1; k2<=(nlstate);k2++){ - for (l2=1; l2<=(nlstate+ndeath);l2++){ - if(l2==k2) continue; - j=(k2-1)*(nlstate+ndeath)+l2; - for (k1=1; k1<=(nlstate);k1++){ - for (l1=1; l1<=(nlstate+ndeath);l1++){ - if(l1==k1) continue; - i=(k1-1)*(nlstate+ndeath)+l1; - if(i<=j) continue; - for (age=bage; age<=fage; age ++){ - if ((int)age %5==0){ - v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM; - v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM; - cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM; - mu1=mu[i][(int) age]/stepm*YEARM ; - mu2=mu[j][(int) age]/stepm*YEARM; - c12=cv12/sqrt(v1*v2); - /* Computing eigen value of matrix of covariance */ - lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; - lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; - if ((lc2 <0) || (lc1 <0) ){ - if(first2==1){ - first1=0; - printf("Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor); - } - fprintf(ficlog,"Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog); - /* lc1=fabs(lc1); */ /* If we want to have them positive */ - /* lc2=fabs(lc2); */ - } + for (i=1; i<=(nlstate)*(nlstate+ndeath);i++) + fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age])); + for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ + fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]); + fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]); + } + i=0; + for (k=1; k<=(nlstate);k++){ + for (l=1; l<=(nlstate+ndeath);l++){ + i++; + fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l); + fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l); + for (j=1; j<=i;j++){ + /* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */ + fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]); + fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age])); + } + } + }/* end of loop for state */ + } /* end of loop for age */ + free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); + free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); + free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); + free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); + + /* Confidence intervalle of pij */ + /* + fprintf(ficgp,"\nunset parametric;unset label"); + fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\""); + fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); + fprintf(fichtm,"\n
    Probability with confidence intervals expressed in year-1 :pijgr%s.png, ",optionfilefiname,optionfilefiname); + fprintf(fichtm,"\n
    ",optionfilefiname); + fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname); + fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob); + */ + + /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ + first1=1;first2=2; + for (k2=1; k2<=(nlstate);k2++){ + for (l2=1; l2<=(nlstate+ndeath);l2++){ + if(l2==k2) continue; + j=(k2-1)*(nlstate+ndeath)+l2; + for (k1=1; k1<=(nlstate);k1++){ + for (l1=1; l1<=(nlstate+ndeath);l1++){ + if(l1==k1) continue; + i=(k1-1)*(nlstate+ndeath)+l1; + if(i<=j) continue; + for (age=bage; age<=fage; age ++){ + if ((int)age %5==0){ + v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM; + v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM; + cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM; + mu1=mu[i][(int) age]/stepm*YEARM ; + mu2=mu[j][(int) age]/stepm*YEARM; + c12=cv12/sqrt(v1*v2); + /* Computing eigen value of matrix of covariance */ + lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + if ((lc2 <0) || (lc1 <0) ){ + if(first2==1){ + first1=0; + printf("Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor); + } + fprintf(ficlog,"Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog); + /* lc1=fabs(lc1); */ /* If we want to have them positive */ + /* lc2=fabs(lc2); */ + } - /* Eigen vectors */ - v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); - /*v21=sqrt(1.-v11*v11); *//* error */ - v21=(lc1-v1)/cv12*v11; - v12=-v21; - v22=v11; - tnalp=v21/v11; - if(first1==1){ - first1=0; - printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); - } - fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); - /*printf(fignu*/ - /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */ - /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */ - if(first==1){ - first=0; - fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n"); - fprintf(ficgp,"\nset parametric;unset label"); - fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); - fprintf(ficgp,"\nset ter svg size 640, 480"); - fprintf(fichtmcov,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ + /* Eigen vectors */ + v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); + /*v21=sqrt(1.-v11*v11); *//* error */ + v21=(lc1-v1)/cv12*v11; + v12=-v21; + v22=v11; + tnalp=v21/v11; + if(first1==1){ + first1=0; + printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); + } + fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); + /*printf(fignu*/ + /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */ + /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */ + if(first==1){ + first=0; + fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n"); + fprintf(ficgp,"\nset parametric;unset label"); + fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); + fprintf(ficgp,"\nset ter svg size 640, 480"); + fprintf(fichtmcov,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ : \ %s_%d%1d%1d-%1d%1d.svg, ",k1,l1,k2,l2,\ - subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \ - subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); - fprintf(fichtmcov,"\n
    ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); - fprintf(fichtmcov,"\n
    Correlation at age %d (%.3f),",(int) age, c12); - fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); - fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); - fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); - fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ - mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ - mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); - }else{ - first=0; - fprintf(fichtmcov," %d (%.3f),",(int) age, c12); - fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); - fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); - fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ - mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ - mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); - }/* if first */ - } /* age mod 5 */ - } /* end loop age */ - fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); - first=1; - } /*l12 */ - } /* k12 */ - } /*l1 */ - }/* k1 */ - } /* loop on combination of covariates j1 */ - free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); - free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); - free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); - free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar); - free_vector(xp,1,npar); - fclose(ficresprob); - fclose(ficresprobcov); - fclose(ficresprobcor); - fflush(ficgp); - fflush(fichtmcov); -} + subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \ + subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); + fprintf(fichtmcov,"\n
    ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); + fprintf(fichtmcov,"\n
    Correlation at age %d (%.3f),",(int) age, c12); + fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); + fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); + fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); + fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ + mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ + mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); + }else{ + first=0; + fprintf(fichtmcov," %d (%.3f),",(int) age, c12); + fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); + fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); + fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ + mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ + mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); + }/* if first */ + } /* age mod 5 */ + } /* end loop age */ + fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); + first=1; + } /*l12 */ + } /* k12 */ + } /*l1 */ + }/* k1 */ + } /* loop on combination of covariates j1 */ + free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); + free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); + free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); + free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar); + free_vector(xp,1,npar); + fclose(ficresprob); + fclose(ficresprobcov); + fclose(ficresprobcor); + fflush(ficgp); + fflush(fichtmcov); + } /******************* Printing html file ***********/ @@ -5536,81 +5539,81 @@ void printinghtml(char fileresu[], char %s
    \n", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); } -fprintf(fichtm," \n
    • Graphs

      • "); + fprintf(fichtm," \n

      • Graphs

        • "); - m=pow(2,cptcoveff); - if (cptcovn < 1) {m=1;ncodemax[1]=1;} + m=pow(2,cptcoveff); + if (cptcovn < 1) {m=1;ncodemax[1]=1;} - jj1=0; - for(k1=1; k1<=m;k1++){ + jj1=0; + for(k1=1; k1<=m;k1++){ - /* for(i1=1; i1<=ncodemax[k1];i1++){ */ - jj1++; - if (cptcovn > 0) { - fprintf(fichtm,"

        ************ Results for covariates"); - for (cpt=1; cpt<=cptcoveff;cpt++){ - fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); - printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); - } - fprintf(fichtm," ************\n
        "); - if(invalidvarcomb[k1]){ - fprintf(fichtm,"\n

        Combination (%d) ignored because no cases

        \n",k1); - printf("\nCombination (%d) ignored because no cases \n",k1); - continue; - } - } - /* 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
        \ + /* for(i1=1; i1<=ncodemax[k1];i1++){ */ + jj1++; + if (cptcovn > 0) { + fprintf(fichtm,"
        ************ Results for covariates"); + for (cpt=1; cpt<=cptcoveff;cpt++){ + fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); + printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); + } + fprintf(fichtm," ************\n
        "); + if(invalidvarcomb[k1]){ + fprintf(fichtm,"\n

        Combination (%d) ignored because no cases

        \n",k1); + printf("\nCombination (%d) ignored because no cases \n",k1); + continue; + } + } + /* 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); - /* 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
        \ + /* 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); - /* Quasi-incidences */ - fprintf(fichtm,"
        \n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ + /* 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); - /* 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
        \ + /* 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); - } - /* 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.\ + } + /* 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); - } - /* 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); - } - if(backcast==1){ - /* Period (stable) back prevalence in each health state */ + } + /* Period (stable) 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
        \ + 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); + } + 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); + } } - } - 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
        \ + if(prevfcast==1){ + /* Projection of prevalence up to period (stable) prevalence in each health state */ + for(cpt=1; cpt<=nlstate;cpt++){ + fprintf(fichtm,"
        \n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s%d_%d.svg
        \ ", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1); - } - } + } + } - for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
        - Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d%d.svg
        \ + 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); - } - /* } /\* end i1 *\/ */ - }/* End k1 */ - fprintf(fichtm,"
      "); + } + /* } /\* end i1 *\/ */ + }/* End k1 */ + fprintf(fichtm,"
    "); - fprintf(fichtm,"\ + fprintf(fichtm,"\ \n

  • Result files (second order: variances)

    \n\ - Parameter file with estimated parameters and covariance matrix: %s
    \ - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file if optimization has been done (mle != 0).
    \ @@ -5622,32 +5625,32 @@ variances but at the covariance matrix. covariance matrix of the one-step probabilities. \ See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres); - fprintf(fichtm," - Standard deviation of one-step probabilities: %s
    \n", - subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_")); - fprintf(fichtm,"\ + fprintf(fichtm," - Standard deviation of one-step probabilities: %s
    \n", + subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_")); + fprintf(fichtm,"\ - Variance-covariance of one-step probabilities: %s
    \n", - subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_")); + subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_")); - fprintf(fichtm,"\ + fprintf(fichtm,"\ - Correlation matrix of one-step probabilities: %s
    \n", - subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_")); - fprintf(fichtm,"\ + subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_")); + fprintf(fichtm,"\ - Variances and covariances of health expectancies by age and initial health status (cov(eij,ekl)(estepm=%2d months): \ %s
    \n
    • ", estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"CVE_")); - fprintf(fichtm,"\ + fprintf(fichtm,"\ - (a) Health expectancies by health status at initial age (eij) and standard errors (in parentheses) (b) life expectancies and standard errors (ei.=ei1+ei2+...)(estepm=%2d months): \ %s
    \n", estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); - fprintf(fichtm,"\ + fprintf(fichtm,"\ - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
    \n", - estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_")); - fprintf(fichtm,"\ + estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_")); + fprintf(fichtm,"\ - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
    \n", - estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_")); - fprintf(fichtm,"\ + estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_")); + fprintf(fichtm,"\ - Standard deviation of period (stable) prevalences: %s
    \n",\ - subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); + subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); /* if(popforecast==1) fprintf(fichtm,"\n */ /* - Prevalences forecasting: f%s
    \n */ @@ -5655,26 +5658,26 @@ See page 'Matrix of variance-covariance /*
    ",fileres,fileres,fileres,fileres); */ /* else */ /* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)

    \n",popforecast, stepm, model); */ - fflush(fichtm); - fprintf(fichtm,"
    • Graphs

      • "); + fflush(fichtm); + fprintf(fichtm,"

      • Graphs

        • "); - m=pow(2,cptcoveff); - if (cptcovn < 1) {m=1;ncodemax[1]=1;} + m=pow(2,cptcoveff); + if (cptcovn < 1) {m=1;ncodemax[1]=1;} - jj1=0; - for(k1=1; k1<=m;k1++){ - /* for(i1=1; i1<=ncodemax[k1];i1++){ */ - jj1++; + jj1=0; + for(k1=1; k1<=m;k1++){ + /* for(i1=1; i1<=ncodemax[k1];i1++){ */ + jj1++; if (cptcovn > 0) { fprintf(fichtm,"

        ************ Results for covariates"); for (cpt=1; cpt<=cptcoveff;cpt++) - fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); + fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); fprintf(fichtm," ************\n
        "); - if(invalidvarcomb[k1]){ - fprintf(fichtm,"\n

        Combination (%d) ignored because no cases

        \n",k1); - continue; - } + if(invalidvarcomb[k1]){ + fprintf(fichtm,"\n

        Combination (%d) ignored because no cases

        \n",k1); + continue; + } } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"\n
        - Observed (cross-sectional) and period (incidence based) \ @@ -5687,10 +5690,10 @@ true period expectancies (those weighted 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); - /* } /\* end i1 *\/ */ - }/* End k1 */ - fprintf(fichtm,"
      "); - fflush(fichtm); + /* } /\* end i1 *\/ */ + }/* End k1 */ + fprintf(fichtm,"
    "); + fflush(fichtm); } /******************* Gnuplot file **************/ @@ -6324,152 +6327,157 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) /*************** Moving average **************/ /* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */ -int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){ + int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){ - int i, cpt, cptcod; - int modcovmax =1; - int mobilavrange, mob; - int iage=0; - - double sum=0.; - double age; - double *sumnewp, *sumnewm; - double *agemingood, *agemaxgood; /* Currently identical for all covariates */ - - - /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ - /* a covariate has 2 modalities, should be equal to ncovcombmax *\/ */ - - sumnewp = vector(1,ncovcombmax); - sumnewm = vector(1,ncovcombmax); - agemingood = vector(1,ncovcombmax); - agemaxgood = vector(1,ncovcombmax); - - for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ - sumnewm[cptcod]=0.; - sumnewp[cptcod]=0.; - agemingood[cptcod]=0; - agemaxgood[cptcod]=0; - } - if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */ + int i, cpt, cptcod; + int modcovmax =1; + int mobilavrange, mob; + int iage=0; + + double sum=0.; + double age; + double *sumnewp, *sumnewm; + double *agemingood, *agemaxgood; /* Currently identical for all covariates */ + + + /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ + /* a covariate has 2 modalities, should be equal to ncovcombmax *\/ */ + + sumnewp = vector(1,ncovcombmax); + sumnewm = vector(1,ncovcombmax); + agemingood = vector(1,ncovcombmax); + agemaxgood = vector(1,ncovcombmax); + + for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ + sumnewm[cptcod]=0.; + sumnewp[cptcod]=0.; + agemingood[cptcod]=0; + agemaxgood[cptcod]=0; + } + if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */ - if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ - if(mobilav==1) mobilavrange=5; /* default */ - else mobilavrange=mobilav; - for (age=bage; age<=fage; age++) - for (i=1; i<=nlstate;i++) - for (cptcod=1;cptcod<=ncovcombmax;cptcod++) - mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; - /* We keep the original values on the extreme ages bage, fage and for - fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 - we use a 5 terms etc. until the borders are no more concerned. - */ - for (mob=3;mob <=mobilavrange;mob=mob+2){ - for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ - for (i=1; i<=nlstate;i++){ - for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ - mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; - for (cpt=1;cpt<=(mob-1)/2;cpt++){ - mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; - mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; - } - mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; - } - } - }/* end age */ - }/* end mob */ - }else - return -1; - for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ - /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */ - agemingood[cptcod]=fage-(mob-1)/2; - for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */ - sumnewm[cptcod]=0.; - for (i=1; i<=nlstate;i++){ - sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; - } - if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ - agemingood[cptcod]=age; - }else{ /* bad */ - for (i=1; i<=nlstate;i++){ - mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; - } /* i */ - } /* end bad */ - }/* age */ - sum=0.; - for (i=1; i<=nlstate;i++){ - sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod]; - } - if(fabs(sum - 1.) > 1.e-3) { /* bad */ - printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod); - /* for (i=1; i<=nlstate;i++){ */ - /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ - /* } /\* i *\/ */ - } /* end bad */ - /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ - /* From youngest, finding the oldest wrong */ - agemaxgood[cptcod]=bage+(mob-1)/2; - for (age=bage+(mob-1)/2; age<=fage; age++){ - sumnewm[cptcod]=0.; - for (i=1; i<=nlstate;i++){ - sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; - } - if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ - agemaxgood[cptcod]=age; - }else{ /* bad */ - for (i=1; i<=nlstate;i++){ - mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; - } /* i */ - } /* end bad */ - }/* age */ - sum=0.; - for (i=1; i<=nlstate;i++){ - sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; - } - if(fabs(sum - 1.) > 1.e-3) { /* bad */ - printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod); - /* for (i=1; i<=nlstate;i++){ */ - /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ - /* } /\* i *\/ */ - } /* end bad */ - - for (age=bage; age<=fage; age++){ - printf("%d %d ", cptcod, (int)age); - sumnewp[cptcod]=0.; - sumnewm[cptcod]=0.; - for (i=1; i<=nlstate;i++){ - sumnewp[cptcod]+=probs[(int)age][i][cptcod]; - sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; - /* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */ - } - /* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */ - } - /* printf("\n"); */ - /* } */ - /* brutal averaging */ - for (i=1; i<=nlstate;i++){ - for (age=1; age<=bage; age++){ - mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; - /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ - } - for (age=fage; age<=AGESUP; age++){ - mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; - /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ - } - } /* end i status */ - for (i=nlstate+1; i<=nlstate+ndeath;i++){ - for (age=1; age<=AGESUP; age++){ - /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/ - mobaverage[(int)age][i][cptcod]=0.; - } - } - }/* end cptcod */ - free_vector(sumnewm,1, ncovcombmax); - free_vector(sumnewp,1, ncovcombmax); - free_vector(agemaxgood,1, ncovcombmax); - free_vector(agemingood,1, ncovcombmax); - return 0; -}/* End movingaverage */ + if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ + if(mobilav==1) mobilavrange=5; /* default */ + else mobilavrange=mobilav; + for (age=bage; age<=fage; age++) + for (i=1; i<=nlstate;i++) + for (cptcod=1;cptcod<=ncovcombmax;cptcod++) + mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; + /* We keep the original values on the extreme ages bage, fage and for + fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 + we use a 5 terms etc. until the borders are no more concerned. + */ + for (mob=3;mob <=mobilavrange;mob=mob+2){ + for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ + for (i=1; i<=nlstate;i++){ + for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ + mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; + for (cpt=1;cpt<=(mob-1)/2;cpt++){ + mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; + mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; + } + mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; + } + } + }/* end age */ + }/* end mob */ + }else + return -1; + for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ + /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */ + if(invalidvarcomb[cptcod]){ + printf("\nCombination (%d) ignored because no cases \n",cptcod); + continue; + } + + agemingood[cptcod]=fage-(mob-1)/2; + for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */ + sumnewm[cptcod]=0.; + for (i=1; i<=nlstate;i++){ + sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; + } + if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ + agemingood[cptcod]=age; + }else{ /* bad */ + for (i=1; i<=nlstate;i++){ + mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; + } /* i */ + } /* end bad */ + }/* age */ + sum=0.; + for (i=1; i<=nlstate;i++){ + sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod]; + } + if(fabs(sum - 1.) > 1.e-3) { /* bad */ + printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod); + /* for (i=1; i<=nlstate;i++){ */ + /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ + /* } /\* i *\/ */ + } /* end bad */ + /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ + /* From youngest, finding the oldest wrong */ + agemaxgood[cptcod]=bage+(mob-1)/2; + for (age=bage+(mob-1)/2; age<=fage; age++){ + sumnewm[cptcod]=0.; + for (i=1; i<=nlstate;i++){ + sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; + } + if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ + agemaxgood[cptcod]=age; + }else{ /* bad */ + for (i=1; i<=nlstate;i++){ + mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; + } /* i */ + } /* end bad */ + }/* age */ + sum=0.; + for (i=1; i<=nlstate;i++){ + sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; + } + if(fabs(sum - 1.) > 1.e-3) { /* bad */ + printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod); + /* for (i=1; i<=nlstate;i++){ */ + /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ + /* } /\* i *\/ */ + } /* end bad */ + + for (age=bage; age<=fage; age++){ + printf("%d %d ", cptcod, (int)age); + sumnewp[cptcod]=0.; + sumnewm[cptcod]=0.; + for (i=1; i<=nlstate;i++){ + sumnewp[cptcod]+=probs[(int)age][i][cptcod]; + sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; + /* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */ + } + /* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */ + } + /* printf("\n"); */ + /* } */ + /* brutal averaging */ + for (i=1; i<=nlstate;i++){ + for (age=1; age<=bage; age++){ + mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; + /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ + } + for (age=fage; age<=AGESUP; age++){ + mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; + /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ + } + } /* end i status */ + for (i=nlstate+1; i<=nlstate+ndeath;i++){ + for (age=1; age<=AGESUP; age++){ + /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/ + mobaverage[(int)age][i][cptcod]=0.; + } + } + }/* end cptcod */ + free_vector(sumnewm,1, ncovcombmax); + free_vector(sumnewp,1, ncovcombmax); + free_vector(agemaxgood,1, ncovcombmax); + free_vector(agemingood,1, ncovcombmax); + return 0; + }/* End movingaverage */ /************** Forecasting ******************/ @@ -8131,6 +8139,10 @@ int hPijx(double *p, int bage, int fage) for(j=1;j<=cptcoveff;j++) fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficrespijb,"******\n"); + if(invalidvarcomb[k]){ + fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k); + continue; + } /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */ for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */ @@ -9648,8 +9660,8 @@ Please run with mle=-1 to get a correct back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); fclose(ficresplb); - /* hBijx(p, bage, fage, mobaverage); */ - /* fclose(ficrespijb); */ + hBijx(p, bage, fage, mobaverage); + fclose(ficrespijb); free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,