/* $Id$
$State$
$Log$
+ 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 ***
/* 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 ***************/
/************* 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
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++){
/* 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++) { */
}
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;
}
}
/************ 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;
+ 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;}
+ /*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;
+ 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<wav[i];mi++){
- m=mw[mi][i];
- agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
- /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
- if(m >=firstpass && m <=lastpass){
- 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] <iagemin-AGEMARGE || (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 (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<wav[i];mi++){
+ m=mw[mi][i];
+ agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
+ /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
+ if(m >=firstpass && m <=lastpass){
+ 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] <iagemin-AGEMARGE || (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 */
+ 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 */
+ /* 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 ***************/
{
/* 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 */
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;
}
}
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));
}
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");
/************ 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<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
- fprintf(fichtm,"\n");
- fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> 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.</li>\n",optionfilehtmcov);
- fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
- fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) 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<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
+ fprintf(fichtm,"\n");
+
+ fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> 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.</li>\n",optionfilehtmcov);
+ fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
+ fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
and drawn. It helps understanding how is the covariance between two incidences.\
They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
- fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
+ fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 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. <br>\
and made the appropriate rotation to look at the uncorrelated principal directions.<br>\
To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.<br> \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<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
+ fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
- 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<h3>Combination (%d) ignored because no cases </h3>\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<h3>Combination (%d) ignored because no cases </h3>\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);
+ 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<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
-
fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
- fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
- fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
- */
+ /* 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<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
+ fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",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); */
- }
+ /* 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<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
+ /* 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<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
:<a href=\"%s_%d%1d%1d-%1d%1d.svg\"> \
%s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\
- subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \
- subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
-
fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
- fprintf(fichtmcov,"\n<br> 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<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+ fprintf(fichtmcov,"\n<br> 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 ***********/
<a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
}
-fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
+
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
- 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,"<hr size=\"2\" color=\"#EC5E5E\">************ 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<hr size=\"2\" color=\"#EC5E5E\">");
- if(invalidvarcomb[k1]){
- fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1);
- printf("\nCombination (%d) ignored because no cases \n",k1);
- continue;
- }
- }
- /* aij, bij */
-
fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
+ /* for(i1=1; i1<=ncodemax[k1];i1++){ */
+ jj1++;
+ if (cptcovn > 0) {
+ fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ 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<hr size=\"2\" color=\"#EC5E5E\">");
+ if(invalidvarcomb[k1]){
+ fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1);
+ printf("\nCombination (%d) ignored because no cases \n",k1);
+ continue;
+ }
+ }
+ /* aij, bij */
+ fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
- /* Pij */
-
fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+ /* Pij */
+ fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
- /* Quasi-incidences */
- fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+ /* Quasi-incidences */
+ fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \
incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \
divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
- /* Survival functions (period) in state j */
- for(cpt=1; cpt<=nlstate;cpt++){
-
fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+ /* Survival functions (period) in state j */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"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,"<br>\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,"<br>\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. \
<a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%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,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
- }
- 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,"<br>\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. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
+ fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
+ }
+ if(backcast==1){
+ /* Period (stable) back prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\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. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
<img src=\"%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,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+ if(prevfcast==1){
+ /* Projection of prevalence up to period (stable) prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
- }
- }
+ }
+ }
- for(cpt=1; cpt<=nlstate;cpt++) {
-
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
+ for(cpt=1; cpt<=nlstate;cpt++) {
+ fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
<img src=\"%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,"</ul>");
+ }
+ /* } /\* end i1 *\/ */
+ }/* End k1 */
+ fprintf(fichtm,"</ul>");
- fprintf(fichtm,"\
+ fprintf(fichtm,"\
\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \
- 95%% confidence intervals and Wald tests of the estimated parameters are in the log file if optimization has been done (mle != 0).<br> \
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: <a href=\"%s\">%s</a> <br>\n",
- subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_"));
- fprintf(fichtm,"\
+
fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
+ subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_"));
+ fprintf(fichtm,"\
- Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
- subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_"));
+ subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_"));
- fprintf(fichtm,"\
+ fprintf(fichtm,"\
- Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\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 <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \
<a href=\"%s\">%s</a> <br>\n</li>",
estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"CVE_"));
- fprintf(fichtm,"\
+ fprintf(fichtm,"\
- (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \
<a href=\"%s\">%s</a> <br>\n</li>",
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), e<sup>ij</sup> 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): <a href=\"%s\">%s</a><br>\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<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\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: <a href=\"%s\">%s</a> <br>\n",\
- subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
+ subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
/* if(popforecast==1) fprintf(fichtm,"\n */
/* - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
/* <br>",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 .)<br><br></li>\n",popforecast, stepm, model); */
- fflush(fichtm);
- fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
+ fflush(fichtm);
+
fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
- 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,"<hr size=\"2\" color=\"#EC5E5E\">************ 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<hr size=\"2\" color=\"#EC5E5E\">");
- if(invalidvarcomb[k1]){
- fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1);
- continue;
- }
+ if(invalidvarcomb[k1]){
+ fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1);
+ continue;
+ }
}
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
drawn in addition to the population based expectancies computed using\
observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\
<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
- /* } /\* end i1 *\/ */
- }/* End k1 */
- fprintf(fichtm,"</ul>");
- fflush(fichtm);
+ /* } /\* end i1 *\/ */
+ }/* End k1 */
+ fprintf(fichtm,"</ul>");
+ fflush(fichtm);
}
/******************* Gnuplot file **************/
/*************** 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 */
+ 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 */
+ /* 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 */
+ 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 */
+ 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 ******************/
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) */
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,
git://henry.ined.fr / .git / commitdiff