********** Variable ");
+ fprintf(ficresphtmfr, "\n
********** Variable ");
+ fprintf(ficlog, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++){
+ if(!FixedV[Tvaraff[z1]]){
+ printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtm, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtmfr, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficlog, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ }else{
+ printf( "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ }
+ }
+ printf( "**********\n#");
fprintf(ficresp, "**********\n#");
+ fprintf(ficresphtm, "**********
\n");
+ fprintf(ficresphtmfr, "**********
\n");
+ fprintf(ficlog, "**********\n");
+ }
+ /*
+ Printing means of quantitative variables if any
+ */
+ for (z1=1; z1<= nqfveff; z1++) {
+ fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.3g (weighted) individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]);
+ fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]);
+ if(weightopt==1){
+ printf(" Weighted mean and standard deviation of");
+ fprintf(ficlog," Weighted mean and standard deviation of");
+ fprintf(ficresphtmfr," Weighted mean and standard deviation of");
+ }
+ /* mu = \frac{w x}{\sum w}
+ var = \frac{\sum w (x-mu)^2}{\sum w} = \frac{w x^2}{\sum w} - mu^2
+ */
+ printf(" fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1]));
+ fprintf(ficlog," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1]));
+ fprintf(ficresphtmfr," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); + } + /* for (z1=1; z1<= nqtveff; z1++) { */ + /* for(m=1;m<=lastpass;m++){ */ + /* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f
\n", z1, m, meanqt[m][z1]); */ + /* } */ + /* } */ + + fprintf(ficresphtm,"
| Age | Prev(%d) | N(%d) | N | ",i,i); } - for(i=1; i<=nlstate;i++) - fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); - fprintf(ficresp, "\n"); + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); + fprintf(ficresphtm, "\n"); - for(i=iagemin; i <= iagemax+3; i++){ - if(i==iagemax+3){ + /* Header of frequency table by age */ + fprintf(ficresphtmfr,"
|---|
| Age | "); + for(s2=-1; s2 <=nlstate+ndeath; s2++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(s2!=0 && m!=0) + fprintf(ficresphtmfr,"%d%d | ",s2,m); + } + } + fprintf(ficresphtmfr, "\n"); + + /* For each age */ + for(iage=iagemin; iage <= iagemax+3; iage++){ + fprintf(ficresphtm,"||||||||
|---|---|---|---|---|---|---|---|---|---|
| 0 | "); + }else if(iage==iagemax+2){ + fprintf(ficlog,"0"); + fprintf(ficresphtmfr,"|||||||||
| Unknown | "); + }else if(iage==iagemax+3){ fprintf(ficlog,"Total"); + fprintf(ficresphtmfr,"|||||||||
| Total | "); }else{ if(first==1){ first=0; printf("See log file for details...\n"); } - fprintf(ficlog,"Age %d", i); + fprintf(ficresphtmfr,"|||||||||
| %d | ",iage); + fprintf(ficlog,"Age %d", iage); } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][i]; + for(s1=1; s1 <=nlstate ; s1++){ + for(m=-1, pp[s1]=0; m <=nlstate+ndeath ; m++) + pp[s1] += freq[s1][m][iage]; } - for(jk=1; jk <=nlstate ; jk++){ + for(s1=1; s1 <=nlstate ; s1++){ for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][i]; - if(pp[jk]>=1.e-10){ + pos += freq[s1][m][iage]; + if(pp[s1]>=1.e-10){ if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + printf(" %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); } - fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); }else{ if(first==1) - printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); + printf(" %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); + fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); } } - - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1,pos=0,posprop=0; jk <=nlstate ; jk++){ - pos += pp[jk]; - posprop += prop[jk][i]; + + for(s1=1; s1 <=nlstate ; s1++){ + /* posprop[s1]=0; */ + for(m=0, pp[s1]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ + pp[s1] += freq[s1][m][iage]; + } /* pp[s1] is the total number of transitions starting from state s1 and any ending status until this age */ + + for(s1=1,pos=0, pospropta=0.; s1 <=nlstate ; s1++){ + pos += pp[s1]; /* pos is the total number of transitions until this age */ + posprop[s1] += prop[s1][iage]; /* prop is the number of transitions from a live state + from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pospropta += prop[s1][iage]; /* prop is the number of transitions from a live state + from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + } + + /* Writing ficresp */ + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + if( iage <= iagemax){ + fprintf(ficresp," %d",iage); + } + }else if( nj==2){ + if( iage <= iagemax){ + fprintf(ficresp," %d",iage); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + } } - for(jk=1; jk <=nlstate ; jk++){ + for(s1=1; s1 <=nlstate ; s1++){ if(pos>=1.e-5){ if(first==1) - printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + printf(" %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos); + fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos); }else{ if(first==1) - printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); + printf(" %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); + fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); } - if( i <= iagemax){ + if( iage <= iagemax){ if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); - /*probs[i][jk][j1]= pp[jk]/pos;*/ - /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + }else if( nj==2){ + fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + } + fprintf(ficresphtm,"%d | %.5f | %.0f | %.0f | ",iage,prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + /*probs[iage][s1][j1]= pp[s1]/pos;*/ + /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/ + } else{ + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); + fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f | ",iage, prop[s1][iage],pospropta); } - else - fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); } - } - - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) - if(freq[jk][m][i] !=0 ) { - if(first==1) - printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]); + pospropt[s1] +=posprop[s1]; + } /* end loop s1 */ + /* pospropt=0.; */ + for(s1=-1; s1 <=nlstate+ndeath; s1++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[s1][m][iage] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); + } + /* printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); */ + fprintf(ficlog," %d%d=%.0f",s1,m,freq[s1][m][iage]); } - if(i <= iagemax) - fprintf(ficresp,"\n"); + if(s1!=0 && m!=0) + fprintf(ficresphtmfr,"%.0f | ",freq[s1][m][iage]); + } + } /* end loop s1 */ + posproptt=0.; + for(s1=1; s1 <=nlstate; s1++){ + posproptt += pospropt[s1]; + } + fprintf(ficresphtmfr,"
| Tot | "); + for(s1=1; s1 <=nlstate ; s1++){ + if(posproptt < 1.e-5){ + fprintf(ficresphtm,"Nanq | %.0f | %.0f | ",pospropt[s1],posproptt); + }else{ + fprintf(ficresphtm,"%.5f | %.0f | %.0f | ",pospropt[s1]/posproptt,pospropt[s1],posproptt); + } + } + fprintf(ficresphtm,"
This combination (%d) is not valid and no result will be produced
",j1); + fprintf(ficresphtmfr,"\nThis combination (%d) is not valid and no result will be produced
",j1); + fprintf(ficlog,"# This combination (%d) is not valid and no result will be produced\n",j1); + printf("# This combination (%d) is not valid and no result will be produced\n",j1); + invalidvarcomb[j1]=1; + }else{ + fprintf(ficresphtm,"\nThis combination (%d) is valid and result will be produced.
",j1); + invalidvarcomb[j1]=0; + } + fprintf(ficresphtmfr,"\n"); + fprintf(ficlog,"\n"); + if(j!=0){ + printf("#Freqsummary: Starting values for combination j1=%d:\n", j1); + for(i=1,s1=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + for(jj=1; jj <=ncovmodel; jj++){ /* For counting s1 */ + if(jj==1){ /* Constant case (in fact cste + age) */ + if(j1==1){ /* All dummy covariates to zero */ + freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */ + freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */ + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3])); + fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); + } + }else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */ + for(iage=iagemin; iage <= iagemax+3; iage++){ + x[iage]= (double)iage; + y[iage]= log(freq[i][k][iage]/freq[i][i][iage]); + /* printf("i=%d, k=%d, s1=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,s1,j1,jj, iage, y[iage]); */ + } + /* Some are not finite, but linreg will ignore these ages */ + no=0; + linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */ + pstart[s1]=b; + pstart[s1-1]=a; + }else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ + printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); + printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); + pstart[s1]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])); + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + printf("s1=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",s1,i,k,s1,p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4])); + }else{ /* Other cases, like quantitative fixed or varying covariates */ + ; + } + /* printf("%12.7f )", param[i][jj][k]); */ + /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ + s1++; + } /* end jj */ + } /* end k!= i */ + } /* end k */ + } /* end i, s1 */ + } /* end j !=0 */ + } /* end selected combination of covariate j1 */ + if(j==0){ /* We can estimate starting values from the occurences in each case */ + printf("#Freqsummary: Starting values for the constants:\n"); + fprintf(ficlog,"\n"); + for(i=1,s1=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + for(jj=1; jj <=ncovmodel; jj++){ + pstart[s1]=p[s1]; /* Setting pstart to p values by default */ + if(jj==1){ /* Age has to be done */ + pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); + printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + } + /* printf("%12.7f )", param[i][jj][k]); */ + /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ + s1++; + } + printf("\n"); + fprintf(ficlog,"\n"); + } + } + } /* end of state i */ + printf("#Freqsummary\n"); + fprintf(ficlog,"\n"); + for(s1=-1; s1 <=nlstate+ndeath; s1++){ + for(s2=-1; s2 <=nlstate+ndeath; s2++){ + /* param[i]|j][k]= freq[s1][s2][iagemax+3] */ + printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); + fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); + /* if(freq[s1][s2][iage] !=0 ) { /\* minimizing output *\/ */ + /* printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ + /* fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ + /* } */ + } + } /* end loop s1 */ + + printf("\n"); + fprintf(ficlog,"\n"); + } /* end j=0 */ + } /* end j */ + + if(mle == -2){ /* We want to use these values as starting values */ + for(i=1, jk=1; i <=nlstate; i++){ + for(j=1; j <=nlstate+ndeath; j++){ + if(j!=i){ + /*ca[0]= k+'a'-1;ca[1]='\0';*/ + printf("%1d%1d",i,j); + fprintf(ficparo,"%1d%1d",i,j); + for(k=1; k<=ncovmodel;k++){ + /* printf(" %lf",param[i][j][k]); */ + /* fprintf(ficparo," %lf",param[i][j][k]); */ + p[jk]=pstart[jk]; + printf(" %f ",pstart[jk]); + fprintf(ficparo," %f ",pstart[jk]); + jk++; + } + printf("\n"); + fprintf(ficparo,"\n"); + } } } - } + } /* end mle=-2 */ dateintmean=dateintsum/k2cpt; - + date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); + fclose(ficresp); - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3); + fclose(ficresphtm); + fclose(ficresphtmfr); + free_vector(idq,1,nqfveff); + free_vector(meanq,1,nqfveff); + free_vector(stdq,1,nqfveff); + free_matrix(meanqt,1,lastpass,1,nqtveff); + free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); + free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); + free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE); + free_vector(pospropt,1,nlstate); + free_vector(posprop,1,nlstate); + free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_vector(pp,1,nlstate); - free_matrix(prop,1,nlstate,iagemin, iagemax+3); - /* End of Freq */ + /* End of freqsummary */ +} + +/* Simple linear regression */ +int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) { + + /* y=a+bx regression */ + double sumx = 0.0; /* sum of x */ + double sumx2 = 0.0; /* sum of x**2 */ + double sumxy = 0.0; /* sum of x * y */ + double sumy = 0.0; /* sum of y */ + double sumy2 = 0.0; /* sum of y**2 */ + double sume2 = 0.0; /* sum of square or residuals */ + double yhat; + + double denom=0; + int i; + int ne=*no; + + for ( i=ifi, ne=0;i<=ila;i++) { + if(!isfinite(x[i]) || !isfinite(y[i])){ + /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */ + continue; + } + ne=ne+1; + sumx += x[i]; + sumx2 += x[i]*x[i]; + sumxy += x[i] * y[i]; + sumy += y[i]; + sumy2 += y[i]*y[i]; + denom = (ne * sumx2 - sumx*sumx); + /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */ + } + + denom = (ne * sumx2 - sumx*sumx); + if (denom == 0) { + // vertical, slope m is infinity + *b = INFINITY; + *a = 0; + if (r) *r = 0; + return 1; + } + + *b = (ne * sumxy - sumx * sumy) / denom; + *a = (sumy * sumx2 - sumx * sumxy) / denom; + if (r!=NULL) { + *r = (sumxy - sumx * sumy / ne) / /* compute correlation coeff */ + sqrt((sumx2 - sumx*sumx/ne) * + (sumy2 - sumy*sumy/ne)); + } + *no=ne; + for ( i=ifi, ne=0;i<=ila;i++) { + if(!isfinite(x[i]) || !isfinite(y[i])){ + /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */ + continue; + } + ne=ne+1; + yhat = y[i] - *a -*b* x[i]; + sume2 += yhat * yhat ; + + denom = (ne * sumx2 - sumx*sumx); + /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */ + } + *sb = sqrt(sume2/(double)(ne-2)/(sumx2 - sumx * sumx /(double)ne)); + *sa= *sb * sqrt(sumx2/ne); + + return 0; } /************ Prevalence ********************/ @@ -2021,134 +5393,243 @@ void prevalence(double ***probs, double We still use firstpass and lastpass as another selection. */ - int i, m, jk, k1, i1, j1, bool, z1,z2,j; - double ***freq; /* Frequencies */ - double *pp, **prop; - double pos,posprop; + int i, m, jk, j1, bool, z1,j, iv; + 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,iagemax+3); + prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ j1=0; - j=cptcoveff; + /*j=cptcoveff;*/ if (cptcovn<1) {j=1;ncodemax[1]=1;} - for(k1=1; k1<=j;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ - j1++; - - for (i=1; i<=nlstate; i++) - for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0.0; - - for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + first=0; + 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+4+AGEMARGE; iage++) + prop[i][iage]=0.0; + printf("Prevalence combination of varying and fixed dummies %d\n",j1); + /* fprintf(ficlog," V%d=%d ",Tvaraff[j1],nbcode[Tvaraff[j1]][codtabm(k,j1)]); */ + fprintf(ficlog,"Prevalence combination of varying and fixed dummies %d\n",j1); + + for (i=1; i<=imx; i++) { /* Each individual */ + bool=1; + /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ + for(mi=1; miComputing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)
%s
\n",digitp); + + varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); + pstamp(ficresvij); + fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); + if(popbased==1) + fprintf(ficresvij,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d\n",mobilav); + else + fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n"); + fprintf(ficresvij,"# Age"); + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate;j++) + fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j); + fprintf(ficresvij,"\n"); + + xp=vector(1,npar); + dnewm=matrix(1,nlstate,1,npar); + doldm=matrix(1,nlstate,1,nlstate); + dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar); + doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); + + gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath); + gpp=vector(nlstate+1,nlstate+ndeath); + gmp=vector(nlstate+1,nlstate+ndeath); + trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ + + if(estepm < stepm){ + printf ("Problem %d lower than %d\n",estepm, stepm); + } + else hstepm=estepm; + /* For example we decided to compute the life expectancy with the smallest unit */ + /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. + nhstepm is the number of hstepm from age to agelim + nstepm is the number of stepm from age to agelim. + Look at function hpijx to understand why because of memory size limitations, + we decided (b) to get a life expectancy respecting the most precise curvature of the + survival function given by stepm (the optimization length). Unfortunately it + means that if the survival funtion is printed every two years of age and if + you sum them up and add 1 year (area under the trapezoids) you won't get the same + results. So we changed our mind and took the option of the best precision. + */ + hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ + agelim = AGESUP; + for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ + nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ + nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + gradg=ma3x(0,nhstepm,1,npar,1,nlstate); + gp=matrix(0,nhstepm,1,nlstate); + gm=matrix(0,nhstepm,1,nlstate); + + + for(theta=1; theta <=npar; theta++){ + for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ + xp[i] = x[i] + (i==theta ?delti[theta]:0); + } + /**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and + * returns into prlim . + */ + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres); + + /* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */ + if (popbased==1) { + if(mobilav ==0){ + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][ij]; + }else{ /* mobilav */ + for(i=1; i<=nlstate;i++) + prlim[i][i]=mobaverage[(int)age][i][ij]; + } + } + /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h. + */ + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */ + /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability + * at horizon h in state j including mortality. + */ + for(j=1; j<= nlstate; j++){ + for(h=0; h<=nhstepm; h++){ + for(i=1, gp[h][j]=0.;i<=nlstate;i++) + gp[h][j] += prlim[i][i]*p3mat[i][j][h]; + } + } + /* Next for computing shifted+ probability of death (h=1 means + computed over hstepm matrices product = hstepm*stepm months) + as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 . + */ + for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(i=1,gpp[j]=0.; i<= nlstate; i++) + gpp[j] += prlim[i][i]*p3mat[i][j][1]; + } + + /* Again with minus shift */ + + for(i=1; i<=npar; i++) /* Computes gradient x - delta */ + xp[i] = x[i] - (i==theta ?delti[theta]:0); + + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres); + + if (popbased==1) { + if(mobilav ==0){ + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][ij]; + }else{ /* mobilav */ + for(i=1; i<=nlstate;i++) + prlim[i][i]=mobaverage[(int)age][i][ij]; + } + } + + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); + + for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ + for(h=0; h<=nhstepm; h++){ + for(i=1, gm[h][j]=0.;i<=nlstate;i++) + gm[h][j] += prlim[i][i]*p3mat[i][j][h]; + } + } + /* This for computing probability of death (h=1 means + computed over hstepm matrices product = hstepm*stepm months) + as a weighted average of prlim. + */ + for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(i=1,gmp[j]=0.; i<= nlstate; i++) + gmp[j] += prlim[i][i]*p3mat[i][j][1]; + } + /* end shifting computations */ + + /**< Computing gradient matrix at horizon h + */ + for(j=1; j<= nlstate; j++) /* vareij */ + for(h=0; h<=nhstepm; h++){ + gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; + } + /**< Gradient of overall mortality p.3 (or p.j) + */ + for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */ + gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; + } + + } /* End theta */ + + /* We got the gradient matrix for each theta and state j */ + trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ + + for(h=0; h<=nhstepm; h++) /* veij */ + for(j=1; j<=nlstate;j++) + for(theta=1; theta <=npar; theta++) + trgradg[h][j][theta]=gradg[h][theta][j]; + + for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ + for(theta=1; theta <=npar; theta++) + trgradgp[j][theta]=gradgp[theta][j]; + /**< as well as its transposed matrix + */ + + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ + for(i=1;i<=nlstate;i++) + for(j=1;j<=nlstate;j++) + vareij[i][j][(int)age] =0.; + + /* Computing trgradg by matcov by gradg at age and summing over h + * and k (nhstepm) formula 15 of article + * Lievre-Brouard-Heathcote + */ + + for(h=0;h<=nhstepm;h++){ + for(k=0;k<=nhstepm;k++){ + matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]); + for(i=1;i<=nlstate;i++) + for(j=1;j<=nlstate;j++) + vareij[i][j][(int)age] += doldm[i][j]*hf*hf; + } + } + + /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of + * p.j overall mortality formula 49 but computed directly because + * we compute the grad (wix pijx) instead of grad (pijx),even if + * wix is independent of theta. + */ + matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); + matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); + for(j=nlstate+1;j<=nlstate+ndeath;j++) + for(i=nlstate+1;i<=nlstate+ndeath;i++) + varppt[j][i]=doldmp[j][i]; + /* end ppptj */ + /* x centered again */ + + prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nres); + + if (popbased==1) { + if(mobilav ==0){ + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][ij]; + }else{ /* mobilav */ + for(i=1; i<=nlstate;i++) + prlim[i][i]=mobaverage[(int)age][i][ij]; + } + } + + /* This for computing probability of death (h=1 means + computed over hstepm (estepm) matrices product = hstepm*stepm months) + as a weighted average of prlim. + */ + hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres); + for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(i=1,gmp[j]=0.;i<= nlstate; i++) + gmp[j] += prlim[i][i]*p3mat[i][j][1]; + } + /* end probability of death */ + + fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); + for(j=nlstate+1; j<=(nlstate+ndeath);j++){ + fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j])); + for(i=1; i<=nlstate;i++){ + fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); + } + } + fprintf(ficresprobmorprev,"\n"); + + fprintf(ficresvij,"%.0f ",age ); + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate;j++){ + fprintf(ficresvij," %.4f", vareij[i][j][(int)age]); + } + fprintf(ficresvij,"\n"); + free_matrix(gp,0,nhstepm,1,nlstate); + free_matrix(gm,0,nhstepm,1,nlstate); + free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate); + free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar); + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + } /* End age */ + free_vector(gpp,nlstate+1,nlstate+ndeath); + free_vector(gmp,nlstate+1,nlstate+ndeath); + free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); + free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ + /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */ + fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480"); + /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ + fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); + fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); + /* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ + /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ + /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ + fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev)); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev)); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev)); + fprintf(fichtm,"\n
File (multiple files are possible if covariates are present): %s\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); + fprintf(fichtm,"\n
Probability is computed over estepm=%d months.
![](\"%s%s.svg\")
\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); + /* fprintf(fichtm,"\n
Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year
![](\"varmuptjgr%s%s.svg\")
\n", stepm,YEARM,digitp,digit); + */ + /* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */ + fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); - if(popbased==1){ - if(mobilav!=0) - strcpy(digitp,"-populbased-mobilav-"); - else strcpy(digitp,"-populbased-nomobil-"); - } - else - strcpy(digitp,"-stablbased-"); - - if (mobilav!=0) { - mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ - fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); - printf(" Error in movingaverage mobilav=%d\n",mobilav); - } - } + free_vector(xp,1,npar); + free_matrix(doldm,1,nlstate,1,nlstate); + free_matrix(dnewm,1,nlstate,1,npar); + free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); + free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar); + free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); + /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */ + fclose(ficresprobmorprev); + fflush(ficgp); + fflush(fichtm); + } /* end varevsij */ - strcpy(fileresprobmorprev,"prmorprev"); - sprintf(digit,"%-d",ij); - /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/ - strcat(fileresprobmorprev,digit); /* Tvar to be done */ - strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */ - strcat(fileresprobmorprev,fileres); - if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) { - printf("Problem with resultfile: %s\n", fileresprobmorprev); - fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); - } - printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); - fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); - fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); - fprintf(ficresprobmorprev,"# Age cov=%-d",ij); - for(j=nlstate+1; j<=(nlstate+ndeath);j++){ - fprintf(ficresprobmorprev," p.%-d SE",j); - for(i=1; i<=nlstate;i++) - fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j); - } - fprintf(ficresprobmorprev,"\n"); - fprintf(ficgp,"\n# Routine varevsij"); - fprintf(fichtm,"\n
Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)
%s
\n",digitp); -/* } */ - varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); +/************ Variance of prevlim ******************/ + void varprevlim(char fileresvpl[], FILE *ficresvpl, double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres) +{ + /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ + /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ - fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n"); - fprintf(ficresvij,"# Age"); + double **dnewmpar,**doldm; + int i, j, nhstepm, hstepm; + double *xp; + double *gp, *gm; + double **gradg, **trgradg; + double **mgm, **mgp; + double age,agelim; + int theta; + + pstamp(ficresvpl); + fprintf(ficresvpl,"# Standard deviation of period (forward stable) prevalences \n"); + fprintf(ficresvpl,"# Age "); + if(nresult >=1) + fprintf(ficresvpl," Result# "); for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate;j++) - fprintf(ficresvij," Cov(e%1d, e%1d)",i,j); - fprintf(ficresvij,"\n"); + fprintf(ficresvpl," %1d-%1d",i,i); + fprintf(ficresvpl,"\n"); xp=vector(1,npar); - dnewm=matrix(1,nlstate,1,npar); + dnewmpar=matrix(1,nlstate,1,npar); doldm=matrix(1,nlstate,1,nlstate); - dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar); - doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); - - gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath); - gpp=vector(nlstate+1,nlstate+ndeath); - gmp=vector(nlstate+1,nlstate+ndeath); - trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ - if(estepm < stepm){ - printf ("Problem %d lower than %d\n",estepm, stepm); - } - else hstepm=estepm; - /* For example we decided to compute the life expectancy with the smallest unit */ - /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. - nhstepm is the number of hstepm from age to agelim - nstepm is the number of stepm from age to agelin. - Look at hpijx to understand the reason of that which relies in memory size - and note for a fixed period like k years */ - /* We decided (b) to get a life expectancy respecting the most precise curvature of the - survival function given by stepm (the optimization length). Unfortunately it - means that if the survival funtion is printed every two years of age and if - you sum them up and add 1 year (area under the trapezoids) you won't get the same - results. So we changed our mind and took the option of the best precision. - */ - hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ + hstepm=1*YEARM; /* Every year of age */ + hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */ agelim = AGESUP; for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ - nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - gradg=ma3x(0,nhstepm,1,npar,1,nlstate); - gp=matrix(0,nhstepm,1,nlstate); - gm=matrix(0,nhstepm,1,nlstate); - + nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ + if (stepm >= YEARM) hstepm=1; + nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ + gradg=matrix(1,npar,1,nlstate); + mgp=matrix(1,npar,1,nlstate); + mgm=matrix(1,npar,1,nlstate); + gp=vector(1,nlstate); + gm=vector(1,nlstate); for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ + for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - - if (popbased==1) { - if(mobilav ==0){ - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav */ - for(i=1; i<=nlstate;i++) - prlim[i][i]=mobaverage[(int)age][i][ij]; - } - } - - for(j=1; j<= nlstate; j++){ - for(h=0; h<=nhstepm; h++){ - for(i=1, gp[h][j]=0.;i<=nlstate;i++) - gp[h][j] += prlim[i][i]*p3mat[i][j][h]; - } + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ + /* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ + /* else */ + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); + for(i=1;i<=nlstate;i++){ + gp[i] = prlim[i][i]; + mgp[theta][i] = prlim[i][i]; } - /* This for computing probability of death (h=1 means - computed over hstepm matrices product = hstepm*stepm months) - as a weighted average of prlim. - */ - for(j=nlstate+1;j<=nlstate+ndeath;j++){ - for(i=1,gpp[j]=0.; i<= nlstate; i++) - gpp[j] += prlim[i][i]*p3mat[i][j][1]; - } - /* end probability of death */ - - for(i=1; i<=npar; i++) /* Computes gradient x - delta */ + for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - - if (popbased==1) { - if(mobilav ==0){ - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav */ - for(i=1; i<=nlstate;i++) - prlim[i][i]=mobaverage[(int)age][i][ij]; - } - } - - for(j=1; j<= nlstate; j++){ - for(h=0; h<=nhstepm; h++){ - for(i=1, gm[h][j]=0.;i<=nlstate;i++) - gm[h][j] += prlim[i][i]*p3mat[i][j][h]; - } - } - /* This for computing probability of death (h=1 means - computed over hstepm matrices product = hstepm*stepm months) - as a weighted average of prlim. - */ - for(j=nlstate+1;j<=nlstate+ndeath;j++){ - for(i=1,gmp[j]=0.; i<= nlstate; i++) - gmp[j] += prlim[i][i]*p3mat[i][j][1]; - } - /* end probability of death */ - - for(j=1; j<= nlstate; j++) /* vareij */ - for(h=0; h<=nhstepm; h++){ - gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; - } - - for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */ - gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ + /* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ + /* else */ + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); + for(i=1;i<=nlstate;i++){ + gm[i] = prlim[i][i]; + mgm[theta][i] = prlim[i][i]; } - + for(i=1;i<=nlstate;i++) + gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; + /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */ } /* End theta */ - trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ - - for(h=0; h<=nhstepm; h++) /* veij */ - for(j=1; j<=nlstate;j++) - for(theta=1; theta <=npar; theta++) - trgradg[h][j][theta]=gradg[h][theta][j]; + trgradg =matrix(1,nlstate,1,npar); - for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ + for(j=1; j<=nlstate;j++) for(theta=1; theta <=npar; theta++) - trgradgp[j][theta]=gradgp[theta][j]; - + trgradg[j][theta]=gradg[theta][j]; + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\nmgm mgp %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf(" %d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\n gradg %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf("%d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf("%d %lf ",theta,gradg[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ - hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1;i<=nlstate;i++) - for(j=1;j<=nlstate;j++) - vareij[i][j][(int)age] =0.; - - for(h=0;h<=nhstepm;h++){ - for(k=0;k<=nhstepm;k++){ - matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]); - for(i=1;i<=nlstate;i++) - for(j=1;j<=nlstate;j++) - vareij[i][j][(int)age] += doldm[i][j]*hf*hf; - } - } - - /* pptj */ - matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); - matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); - for(j=nlstate+1;j<=nlstate+ndeath;j++) - for(i=nlstate+1;i<=nlstate+ndeath;i++) - varppt[j][i]=doldmp[j][i]; - /* end ppptj */ - /* x centered again */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); - - if (popbased==1) { - if(mobilav ==0){ - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav */ - for(i=1; i<=nlstate;i++) - prlim[i][i]=mobaverage[(int)age][i][ij]; - } + varpl[i][(int)age] =0.; + if((int)age==79 ||(int)age== 80 ||(int)age== 81){ + matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg); + }else{ + matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg); } - - /* This for computing probability of death (h=1 means - computed over hstepm (estepm) matrices product = hstepm*stepm months) - as a weighted average of prlim. - */ - for(j=nlstate+1;j<=nlstate+ndeath;j++){ - for(i=1,gmp[j]=0.;i<= nlstate; i++) - gmp[j] += prlim[i][i]*p3mat[i][j][1]; - } - /* end probability of death */ - - fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); - for(j=nlstate+1; j<=(nlstate+ndeath);j++){ - fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j])); - for(i=1; i<=nlstate;i++){ - fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); - } - } - fprintf(ficresprobmorprev,"\n"); + for(i=1;i<=nlstate;i++) + varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ - fprintf(ficresvij,"%.0f ",age ); - for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate;j++){ - fprintf(ficresvij," %.4f", vareij[i][j][(int)age]); - } - fprintf(ficresvij,"\n"); - free_matrix(gp,0,nhstepm,1,nlstate); - free_matrix(gm,0,nhstepm,1,nlstate); - free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate); - free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar); - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + fprintf(ficresvpl,"%.0f ",age ); + if(nresult >=1) + fprintf(ficresvpl,"%d ",nres ); + for(i=1; i<=nlstate;i++){ + fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); + /* for(j=1;j<=nlstate;j++) */ + /* fprintf(ficresvpl," %d %.5f ",j,prlim[j][i]); */ + } + fprintf(ficresvpl,"\n"); + free_vector(gp,1,nlstate); + free_vector(gm,1,nlstate); + free_matrix(mgm,1,npar,1,nlstate); + free_matrix(mgp,1,npar,1,nlstate); + free_matrix(gradg,1,npar,1,nlstate); + free_matrix(trgradg,1,nlstate,1,npar); } /* End age */ - free_vector(gpp,nlstate+1,nlstate+ndeath); - free_vector(gmp,nlstate+1,nlstate+ndeath); - free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); - free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ - fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65"); - /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ - fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); -/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ -/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ -/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ - fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l 1 ",subdirf(fileresprobmorprev)); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",subdirf(fileresprobmorprev)); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l 2 ",subdirf(fileresprobmorprev)); - fprintf(fichtm,"\n
File (multiple files are possible if covariates are present): %s\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); - fprintf(fichtm,"\n
Probability is computed over estepm=%d months.
![](\"%s%s.png\")
\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); - /* fprintf(fichtm,"\n
Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year
![](\"varmuptjgr%s%s.png\")
\n", stepm,YEARM,digitp,digit); -*/ -/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit); */ - fprintf(ficgp,"\nset out \"%s%s.png\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); free_vector(xp,1,npar); - free_matrix(doldm,1,nlstate,1,nlstate); - free_matrix(dnewm,1,nlstate,1,npar); - free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); - free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar); - free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); - if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - fclose(ficresprobmorprev); - fflush(ficgp); - fflush(fichtm); -} /* end varevsij */ + free_matrix(doldm,1,nlstate,1,npar); + free_matrix(dnewmpar,1,nlstate,1,nlstate); -/************ Variance of prevlim ******************/ -void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) +} + + +/************ Variance of backprevalence limit ******************/ + void varbrevlim(char fileresvbl[], FILE *ficresvbl, double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres) { - /* Variance of prevalence limit */ + /* Variance of backward prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ - double **newm; - double **dnewm,**doldm; + + double **dnewmpar,**doldm; int i, j, nhstepm, hstepm; - int k, cptcode; double *xp; double *gp, *gm; double **gradg, **trgradg; + double **mgm, **mgp; double age,agelim; int theta; - - fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); - fprintf(ficresvpl,"# Age"); + + pstamp(ficresvbl); + fprintf(ficresvbl,"# Standard deviation of back (stable) prevalences \n"); + fprintf(ficresvbl,"# Age "); + if(nresult >=1) + fprintf(ficresvbl," Result# "); for(i=1; i<=nlstate;i++) - fprintf(ficresvpl," %1d-%1d",i,i); - fprintf(ficresvpl,"\n"); + fprintf(ficresvbl," %1d-%1d",i,i); + fprintf(ficresvbl,"\n"); xp=vector(1,npar); - dnewm=matrix(1,nlstate,1,npar); + dnewmpar=matrix(1,nlstate,1,npar); doldm=matrix(1,nlstate,1,nlstate); hstepm=1*YEARM; /* Every year of age */ hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */ - agelim = AGESUP; - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ + agelim = AGEINF; + for (age=fage; age>=bage; age --){ /* If stepm=6 months */ + nhstepm=(int) rint((age-agelim)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ if (stepm >= YEARM) hstepm=1; nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ gradg=matrix(1,npar,1,nlstate); + mgp=matrix(1,npar,1,nlstate); + mgm=matrix(1,npar,1,nlstate); gp=vector(1,nlstate); gm=vector(1,nlstate); @@ -2790,18 +6739,27 @@ void varprevlim(char fileres[], double * for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - for(i=1;i<=nlstate;i++) - gp[i] = prlim[i][i]; - - for(i=1; i<=npar; i++) /* Computes gradient */ + if(mobilavproj > 0 ) + bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres); + else + bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres); + for(i=1;i<=nlstate;i++){ + gp[i] = bprlim[i][i]; + mgp[theta][i] = bprlim[i][i]; + } + for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - for(i=1;i<=nlstate;i++) - gm[i] = prlim[i][i]; - + if(mobilavproj > 0 ) + bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres); + else + bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres); + for(i=1;i<=nlstate;i++){ + gm[i] = bprlim[i][i]; + mgm[theta][i] = bprlim[i][i]; + } for(i=1;i<=nlstate;i++) gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; + /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */ } /* End theta */ trgradg =matrix(1,nlstate,1,npar); @@ -2809,112 +6767,140 @@ void varprevlim(char fileres[], double * for(j=1; j<=nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[j][theta]=gradg[theta][j]; + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\nmgm mgp %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf(" %d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\n gradg %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf("%d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf("%d %lf ",theta,gradg[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ for(i=1;i<=nlstate;i++) - varpl[i][(int)age] =0.; - matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + varbpl[i][(int)age] =0.; + if((int)age==79 ||(int)age== 80 ||(int)age== 81){ + matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg); + }else{ + matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg); + } for(i=1;i<=nlstate;i++) - varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ + varbpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ - fprintf(ficresvpl,"%.0f ",age ); + fprintf(ficresvbl,"%.0f ",age ); + if(nresult >=1) + fprintf(ficresvbl,"%d ",nres ); for(i=1; i<=nlstate;i++) - fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); - fprintf(ficresvpl,"\n"); + fprintf(ficresvbl," %.5f (%.5f)",bprlim[i][i],sqrt(varbpl[i][(int)age])); + fprintf(ficresvbl,"\n"); free_vector(gp,1,nlstate); free_vector(gm,1,nlstate); + free_matrix(mgm,1,npar,1,nlstate); + free_matrix(mgp,1,npar,1,nlstate); free_matrix(gradg,1,npar,1,nlstate); free_matrix(trgradg,1,nlstate,1,npar); } /* End age */ free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,npar); - free_matrix(dnewm,1,nlstate,1,nlstate); + free_matrix(dnewmpar,1,nlstate,1,nlstate); } /************ 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) -{ - int i, j=0, i1, k1, l1, t, tj; - int k2, l2, j1, z1; - int k=0,l, cptcode; - int first=1, first1; - 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,agelim, cov[NCOVMAX]; - 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,fileres); - 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,fileres); - 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); - - fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); - fprintf(ficresprob,"# Age"); - fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); - fprintf(ficresprobcov,"# Age"); - fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); - fprintf(ficresprobcov,"# 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
Matrix of variance-covariance of pairs of step probabilities (drawings)
Matrix of variance-covariance of pairs of step probabilities
\n\ - file %s\n",optionfilehtmcov); - fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p
Computing and drawing one step probabilities with their confidence intervals
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. File %s\n\n
Matrix of variance-covariance of pairs of step probabilities
\n",optionfilehtmcov, optionfilehtmcov); + fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p\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.
\ @@ -2922,310 +6908,521 @@ 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; - if (cptcovn<1) {tj=1;ncodemax[1]=1;} - j1=0; - for(t=1; t<=tj;t++){ - for(i1=1; i1<=ncodemax[t];i1++){ - j1++; - if (cptcovn>0) { - fprintf(ficresprob, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[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]][codtab[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]][codtab[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]][codtab[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]][codtab[j1][z1]]); - fprintf(ficresprobcor, "**********\n#"); - } - - for (age=bage; age<=fage; age ++){ - cov[2]=age; - for (k=1; k<=cptcovn;k++) { - cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]]; - } - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; - for (k=1; k<=cptcovprod;k++) - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; - - 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(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); - - 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); - - 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(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); - 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); - - 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]; - - /*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); - - 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=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++){ - 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 */ - - /* Confidence intervalle of pij */ - /* - fprintf(ficgp,"\nset noparametric;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
![](\"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);
- */
+ 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 or only once*/
+ for(nres=1;nres <=1; nres++){ /* For each resultline */
+ /* for(nres=1;nres <=nresult; nres++){ /\* For each resultline *\/ */
+ 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(fichtmcov, "\n********** Variable "); + /* for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); */ + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "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)]; */ + for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ + /* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[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;*/ + } + /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */ + /* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] */ + /*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + for (k=1; k<=cptcovage;k++){ /* For product with age */ + if(Dummy[Tage[k]]==2){ /* dummy with age */ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,k)]*cov[2]; + /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ + } else if(Dummy[Tage[k]]==3){ /* quantitative with age */ + printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]); + exit(1); + /* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\* Using the mean of quantitative variable Tvar[Tage[k]] /\* Tqresult[nres][k]; *\/ */ + /* cov[++k1]=Tqresult[nres][k]; */ + } + /* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ + } + for (k=1; k<=cptcovprod;k++){/* For product without age */ + if(Dummy[Tvard[k][1]==0]){ + if(Dummy[Tvard[k][2]==0]){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,k)] * nbcode[Tvard[k][2]][codtabm(j1,k)]; + /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ + }else{ /* Should we use the mean of the quantitative variables? */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,k)] * Tqresult[nres][k]; + /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */ + } + }else{ + if(Dummy[Tvard[k][2]==0]){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(j1,k)] * Tqinvresult[nres][Tvard[k][1]]; + /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */ + }else{ + cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; + /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; */ + } + } + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */ + } +/* For each age and combination of dummy covariates we slightly move the parameters of delti in order to get the gradient*/ + 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); + + 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); + + 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]; + } - /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ - first1=1; - 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.; - /* 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,"\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 png small\nset size 0.65,0.65"); - fprintf(fichtmcov,"\nEllipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ - :\ -%s%d%1d%1d-%1d%1d.png, ",k1,l1,k2,l2,\ - subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\ - subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - fprintf(fichtmcov,"\n
![](\"%s%d%1d%1d-%1d%1d.png\")
",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
- fprintf(fichtmcov,"\nCorrelation at age %d (%.3f),",(int) age, c12); - fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",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 \"%s%d%1d%1d-%1d%1d.png\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - first=1; - } /*l12 */ - } /* k12 */ - } /*l1 */ - }/* k1 */ - } /* loop covariates */ - } - 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_vector(xp,1,npar); - fclose(ficresprob); - fclose(ficresprobcov); - fclose(ficresprobcor); - fflush(ficgp); - fflush(fichtmcov); -} + 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); + + 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]; + + /*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); + + 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 */
+ if(1+(v1-lc1)*(v1-lc1)/cv12/cv12 <1.e-5){
+ printf(" Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
+ fprintf(ficlog," Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
+ v11=(1./sqrt(fabs(1+(v1-lc1)*(v1-lc1)/cv12/cv12)));
+ }else
+ 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![](\"%s_%d%1d%1d-%1d%1d.svg\")
",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(fabs(lc1)),v12,sqrt(fabs(lc2)), \
+ mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */
+ }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(fabs(lc2)), \
+ mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(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 nres */
+ } /* 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 ***********/
-void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
+void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
int lastpass, int stepm, int weightopt, char model[],\
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
- int popforecast, int estepm ,\
- double jprev1, double mprev1,double anprev1, \
- double jprev2, double mprev2,double anprev2){
- int jj1, k1, i1, cpt;
-
- fprintf(fichtm,"
Result files (first order: no variance)
\n \ - - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s
\n ", - jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p")); + int popforecast, int mobilav, int prevfcast, int mobilavproj, int prevbcast, int estepm , \ + double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \ + double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){ + int jj1, k1, i1, cpt, k4, nres; + /* In fact some results are already printed in fichtm which is open */ + fprintf(fichtm,""); +/* fprintf(fichtm,"- model=1+age+%s\n \ */ +/*
Result files (first order: no variance)
\n"); + fprintf(fichtm,"- - Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s (html file)
\n", + jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm")); + fprintf(fichtm," - - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s (html file) ",
+ jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm"));
+ fprintf(fichtm,", %s (text file)
\n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_")); fprintf(fichtm,"\ - Estimated transition probabilities over %d (stepm) months: %s
\n ", - stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij")); + stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_")); fprintf(fichtm,"\ - - Stable prevalence in each health state: %s
\n", - subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); + - Estimated back transition probabilities over %d (stepm) months: %s
\n ", + stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_")); fprintf(fichtm,"\ - - Life expectancies by age and initial health status (estepm=%2d months): \ - %s
\n ",
- estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
+ - Period (forward) prevalence in each health state: %s
\n", + subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); + fprintf(fichtm,"\ + - Backward prevalence in each health state: %s
\n", + subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_")); + fprintf(fichtm,"\ + - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ + %s
\n", + estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); + if(prevfcast==1){ + fprintf(fichtm,"\ + - Prevalence projections by age and states: \ + %s
\n ", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
+ }
+
+
+ m=pow(2,cptcoveff);
+ if (cptcovn < 1) {m=1;ncodemax[1]=1;}
+
+ fprintf(fichtm," \n- Graphs (first order)
- "); + + /* if(nqfveff+nqtveff 0) */ /* Test to be done */ + fprintf(fichtm,"************ Results for covariates"); + for (cpt=1; cpt<=cptcoveff;cpt++){ + fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]); + } + for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ + fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + } + if(invalidvarcomb[k1]){ + fprintf(fichtm," Warning Combination (%d) ignored because no cases ",k1); + continue; + } + fprintf(fichtm," "); + } /* cptcovn >0 */ + } + fprintf(fichtm," \n
- Graphs
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).
\ +But because parameters are usually highly correlated (a higher incidence of disability \ +and a higher incidence of recovery can give very close observed transition) it might \ +be very useful to look not only at linear confidence intervals estimated from the \ +variances but at the covariance matrix. And instead of looking at the estimated coefficients \ +(parameters) of the logistic regression, it might be more meaningful to visualize the \ +covariance matrix of the one-step probabilities. \ +See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres); - fprintf(fichtm,"\ -\nResult files (second order: variances)
\n\ - - Parameter file with estimated parameters and covariance matrix: %s
\n", rfileres,rfileres); - - fprintf(fichtm," - Variance of one-step probabilities: %s
\n", - subdirf2(fileres,"prob"),subdirf2(fileres,"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(fileres,"probcov"),subdirf2(fileres,"probcov")); + subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_")); - fprintf(fichtm,"\ + fprintf(fichtm,"\ - Correlation matrix of one-step probabilities: %s
\n", - subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor")); - fprintf(fichtm,"\ - - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): %s
\n", - estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); - fprintf(fichtm,"\ - - Health expectancies with their variances (no covariance): %s
\n", - subdirf2(fileres,"t"),subdirf2(fileres,"t")); - fprintf(fichtm,"\ - - Standard deviation of stable prevalences: %s
\n",\ - subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); + 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,"\
+ - (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 - Graphs
/*
",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
- Graphs (second order)
- "); + + /* if(nqfveff+nqtveff 0) */ /* Test to be done */ + fprintf(fichtm,"************ Results for covariates"); + for (cpt=1; cpt<=cptcoveff;cpt++){ + fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]); + } + for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ + fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + } + if(invalidvarcomb[k1]){ + fprintf(fichtm," Warning Combination (%d) ignored because no cases ",k1); + continue; + } + fprintf(fichtm," "); + } /* cptcovn >0 */ + } + fprintf(fichtm," \n
Result files
\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):
"); + fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year
",p[1],p[2],agegomp); + for (i=1;i<=2;i++) + fprintf(fichtm," p[%d] = %lf [%f ; %f]
\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); + fprintf(fichtm,"![](\"graphmort.svg\")
");
+ fprintf(fichtm,"Life table
\n
"); + + fprintf(fichtm,"\nAge lx qx d(x,x+1) Lx Tx e
"); + + for (k=agegomp;k<(agemortsup-2);k++) + fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf
\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]); + + + fflush(fichtm); +} + +/******************* Gnuplot file **************/ +void printinggnuplotmort(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ + + char dirfileres[132],optfileres[132]; + + int ng; + + + /*#ifdef windows */ + fprintf(ficgp,"cd \"%s\" \n",pathc); + /*#endif */ + + + strcpy(dirfileres,optionfilefiname); + strcpy(optfileres,"vpl"); + fprintf(ficgp,"set out \"graphmort.svg\"\n "); + fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); + fprintf(ficgp, "set ter svg size 640, 480\n set log y\n"); + /* fprintf(ficgp, "set size 0.65,0.65\n"); */ + fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp); + +} + +int readdata(char datafile[], int firstobs, int lastobs, int *imax) +{ + + /*-------- data file ----------*/ + FILE *fic; + char dummy[]=" "; + int i=0, j=0, n=0, iv=0, v; + int lstra; + int linei, month, year,iout; + int noffset=0; /* This is the offset if BOM data file */ + char line[MAXLINE], linetmp[MAXLINE]; + char stra[MAXLINE], strb[MAXLINE]; + char *stratrunc; + + DummyV=ivector(1,NCOVMAX); /* 1 to 3 */ + FixedV=ivector(1,NCOVMAX); /* 1 to 3 */ + + for(v=1; v <=ncovcol;v++){ + DummyV[v]=0; + FixedV[v]=0; + } + for(v=ncovcol+1; v <=ncovcol+nqv;v++){ + DummyV[v]=1; + FixedV[v]=0; + } + for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){ + DummyV[v]=0; + FixedV[v]=1; + } + for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){ + DummyV[v]=1; + FixedV[v]=1; + } + for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){ + printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); + fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); + } + + if((fic=fopen(datafile,"r"))==NULL) { + printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout); + fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1; + } + + /* Is it a BOM UTF-8 Windows file? */ + /* First data line */ + linei=0; + while(fgets(line, MAXLINE, fic)) { + noffset=0; + if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */ + { + noffset=noffset+3; + printf("# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile);fflush(stdout); + fprintf(ficlog,"# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile); + fflush(ficlog); return 1; + } + /* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/ + else if( line[0] == (char)0xFF && line[1] == (char)0xFE) + { + noffset=noffset+2; + printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); + fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); + fflush(ficlog); return 1; + } + else if( line[0] == 0 && line[1] == 0) + { + if( line[2] == (char)0xFE && line[3] == (char)0xFF){ + noffset=noffset+4; + printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); + fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); + fflush(ficlog); return 1; + } + } else{ + ;/*printf(" Not a BOM file\n");*/ + } + /* If line starts with a # it is a comment */ + if (line[noffset] == '#') { + linei=linei+1; + break; + }else{ + break; + } + } + fclose(fic); + if((fic=fopen(datafile,"r"))==NULL) { + printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout); + fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1; + } + /* Not a Bom file */ + + i=1; + while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { + linei=linei+1; + for(j=strlen(line); j>=0;j--){ /* Untabifies line */ + if(line[j] == '\t') + line[j] = ' '; + } + for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){ + ; + }; + line[j+1]=0; /* Trims blanks at end of line */ + if(line[0]=='#'){ + fprintf(ficlog,"Comment line\n%s\n",line); + printf("Comment line\n%s\n",line); + continue; + } + trimbb(linetmp,line); /* Trims multiple blanks in line */ + strcpy(line, linetmp); + + /* Loops on waves */ + for (j=maxwav;j>=1;j--){ + for (iv=nqtv;iv>=1;iv--){ /* Loop on time varying quantitative variables */ + cutv(stra, strb, line, ' '); + if(strb[0]=='.') { /* Missing value */ + lval=-1; + cotqvar[j][iv][i]=-1; /* 0.0/0.0 */ + cotvar[j][ntv+iv][i]=-1; /* For performance reasons */ + if(isalpha(strb[1])) { /* .m or .d Really Missing value */ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog); + return 1; + } + }else{ + errno=0; + /* what_kind_of_number(strb); */ + dval=strtod(strb,&endptr); + /* if( strb[0]=='\0' || (*endptr != '\0')){ */ + /* if(strb != endptr && *endptr == '\0') */ + /* dval=dlval; */ + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog); + return 1; + } + cotqvar[j][iv][i]=dval; + cotvar[j][ntv+iv][i]=dval; + } + strcpy(line,stra); + }/* end loop ntqv */ + + for (iv=ntv;iv>=1;iv--){ /* Loop on time varying dummies */ + cutv(stra, strb, line, ' '); + if(strb[0]=='.') { /* Missing value */ + lval=-1; + }else{ + errno=0; + lval=strtol(strb,&endptr,10); + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog); + return 1; + } + } + if(lval <-1 || lval >1){ + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + Should be a value of %d(nth) covariate of wave %d (0 should be the value for the reference and 1\n \ + for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ + output of IMaCh is often meaningless.\n \ + Exiting.\n",lval,linei, i,line,iv,j); + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + Should be a value of %d(nth) covariate of wave %d (0 should be the value for the reference and 1\n \ + for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ + output of IMaCh is often meaningless.\n \ + Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog); + return 1; + } + cotvar[j][iv][i]=(double)(lval); + strcpy(line,stra); + }/* end loop ntv */ + + /* Statuses at wave */ + cutv(stra, strb, line, ' '); + if(strb[0]=='.') { /* Missing value */ + lval=-1; + }else{ + errno=0; + lval=strtol(strb,&endptr,10); + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog); + return 1; + } + } + + s[j][i]=lval; + + /* Date of Interview */ + strcpy(line,stra); + cutv(stra, strb,line,' '); + if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ + } + else if( (iout=sscanf(strb,"%s.",dummy)) != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog); + return 1; + } + anint[j][i]= (double) year; + mint[j][i]= (double)month; + /* if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ */ + /* printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */ + /* fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */ + /* } */ + strcpy(line,stra); + } /* End loop on waves */ + + /* Date of death */ + cutv(stra, strb,line,' '); + if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ + } + else if( (iout=sscanf(strb,"%s.",dummy)) != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + return 1; + } + andc[i]=(double) year; + moisdc[i]=(double) month; + strcpy(line,stra); + + /* Date of birth */ + cutv(stra, strb,line,' '); + if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ + } + else if( (iout=sscanf(strb,"%s.", dummy)) != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + return 1; + } + if (year==9999) { + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog); + return 1; + + } + annais[i]=(double)(year); + moisnais[i]=(double)(month); + for (j=1;j<=maxwav;j++){ + if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ + printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j,(int)moisnais[i],(int)annais[i]); + fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j, (int)moisnais[i],(int)annais[i]); + } + } + + strcpy(line,stra); + + /* Sample weight */ + cutv(stra, strb,line,' '); + errno=0; + dval=strtod(strb,&endptr); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); + fprintf(ficlog,"Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); + fflush(ficlog); + return 1; + } + weight[i]=dval; + strcpy(line,stra); + + for (iv=nqv;iv>=1;iv--){ /* Loop on fixed quantitative variables */ + cutv(stra, strb, line, ' '); + if(strb[0]=='.') { /* Missing value */ + lval=-1; + coqvar[iv][i]=NAN; + covar[ncovcol+iv][i]=NAN; /* including qvar in standard covar for performance reasons */ + }else{ + errno=0; + /* what_kind_of_number(strb); */ + dval=strtod(strb,&endptr); + /* if(strb != endptr && *endptr == '\0') */ + /* dval=dlval; */ + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value (out of %d) constant for all waves. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqv, maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value (out of %d) constant for all waves. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqv, maxwav);fflush(ficlog); + return 1; + } + coqvar[iv][i]=dval; + covar[ncovcol+iv][i]=dval; /* including qvar in standard covar for performance reasons */ + } + strcpy(line,stra); + }/* end loop nqv */ + + /* Covariate values */ + for (j=ncovcol;j>=1;j--){ + cutv(stra, strb,line,' '); + if(strb[0]=='.') { /* Missing covariate value */ + lval=-1; + }else{ + errno=0; + lval=strtol(strb,&endptr,10); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line); + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog); + return 1; + } + } + if(lval <-1 || lval >1){ + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ + for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ + output of IMaCh is often meaningless.\n \ + Exiting.\n",lval,linei, i,line,j); + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ + for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ + output of IMaCh is often meaningless.\n \ + Exiting.\n",lval,linei, i,line,j);fflush(ficlog); + return 1; + } + covar[j][i]=(double)(lval); + strcpy(line,stra); + } + lstra=strlen(stra); + + if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ + stratrunc = &(stra[lstra-9]); + num[i]=atol(stratrunc); + } + else + num[i]=atol(stra); + /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ + printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ + + i=i+1; + } /* End loop reading data */ + + *imax=i-1; /* Number of individuals */ + fclose(fic); + + return (0); + /* endread: */ + printf("Exiting readdata: "); + fclose(fic); + return (1); +} + +void removefirstspace(char **stri){/*, char stro[]) {*/ + char *p1 = *stri, *p2 = *stri; + while (*p2 == ' ') + p2++; + /* while ((*p1++ = *p2++) !=0) */ + /* ; */ + /* do */ + /* while (*p2 == ' ') */ + /* p2++; */ + /* while (*p1++ == *p2++); */ + *stri=p2; +} + +int decoderesult ( char resultline[], int nres) +/**< This routine decode one result line and returns the combination # of dummy covariates only **/ +{ + int j=0, k=0, k1=0, k2=0, k3=0, k4=0, match=0, k2q=0, k3q=0, k4q=0; + char resultsav[MAXLINE]; + int resultmodel[MAXLINE]; + int modelresult[MAXLINE]; + char stra[80], strb[80], strc[80], strd[80],stre[80]; + + removefirstspace(&resultline); + + if (strstr(resultline,"v") !=0){ + printf("Error. 'v' must be in upper case 'V' result: %s ",resultline); + fprintf(ficlog,"Error. 'v' must be in upper case result: %s ",resultline);fflush(ficlog); + return 1; + } + trimbb(resultsav, resultline); + if (strlen(resultsav) >1){ + j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */ + } + if(j == 0){ /* Resultline but no = */ + TKresult[nres]=0; /* Combination for the nresult and the model */ + return (0); + } + if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ + printf("ERROR: the number of variables in this result line, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); + fprintf(ficlog,"ERROR: the number of variables in the resultline, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); + } + for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */ + if(nbocc(resultsav,'=') >1){ + cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//* resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */ + cutl(strc,strd,strb,'='); /* strb:"V4=1" strc="1" strd="V4" */ + }else + cutl(strc,strd,resultsav,'='); + Tvalsel[k]=atof(strc); /* 1 */ /* Tvalsel of k is the float value of the kth covariate appearing in this result line */ + + cutl(strc,stre,strd,'V'); /* strd='V4' strc=4 stre='V' */; + Tvarsel[k]=atoi(strc); /* 4 */ /* Tvarsel is the id of the kth covariate in the result line Tvarsel[1] in "V4=1.." is 4.*/ + /* Typevarsel[k]=1; /\* 1 for age product *\/ */ + /* cptcovsel++; */ + if (nbocc(stra,'=') >0) + strcpy(resultsav,stra); /* and analyzes it */ + } + /* Checking for missing or useless values in comparison of current model needs */ + for(k1=1; k1<= cptcovt ;k1++){ /* Loop on model. model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + if(Typevar[k1]==0){ /* Single covariate in model *//*0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ + match=0; + for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ + if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5 */ + modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2 modelresult[3]=3 modelresult[6]=4 modelresult[9]=5 */ + match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */ + break; + } + } + if(match == 0){ + printf("Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); + fprintf(ficlog,"Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); + return 1; + } + } + } + /* Checking for missing or useless values in comparison of current model needs */ + for(k2=1; k2 <=j;k2++){ /* Loop on resultline variables: result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ + match=0; + for(k1=1; k1<= cptcovt ;k1++){ /* loop on model: model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + if(Typevar[k1]==0){ /* Single */ + if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4 */ + resultmodel[k1]=k2; /* k2th variable of the model corresponds to k1 variable of the model. resultmodel[2]=1 resultmodel[1]=2 resultmodel[3]=3 resultmodel[6]=4 resultmodel[9]=5 */ + ++match; + } + } + } + if(match == 0){ + printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + fprintf(ficlog,"Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + return 1; + }else if(match > 1){ + printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); + fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); + return 1; + } + } + + /* We need to deduce which combination number is chosen and save quantitative values */ + /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + /* result line V4=1 V5=25.1 V3=0 V2=8 V1=1 */ + /* should give a combination of dummy V4=1, V3=0, V1=1 => V4*2**(0) + V3*2**(1) + V1*2**(2) = 5 + (1offset) = 6*/ + /* result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ + /* should give a combination of dummy V4=1, V3=1, V1=0 => V4*2**(0) + V3*2**(1) + V1*2**(2) = 3 + (1offset) = 4*/ + /* 1 0 0 0 */ + /* 2 1 0 0 */ + /* 3 0 1 0 */ + /* 4 1 1 0 */ /* V4=1, V3=1, V1=0 */ + /* 5 0 0 1 */ + /* 6 1 0 1 */ /* V4=1, V3=0, V1=1 */ + /* 7 0 1 1 */ + /* 8 1 1 1 */ + /* V(Tvresult)=Tresult V4=1 V3=0 V1=1 Tresult[nres=1][2]=0 */ + /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */ + /* V5*age V5 known which value for nres? */ + /* Tqinvresult[2]=8 Tqinvresult[1]=25.1 */ + for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* loop on model line */ + if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */ + k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */ + k2=(int)Tvarsel[k3]; /* Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */ + k+=Tvalsel[k3]*pow(2,k4); /* Tvalsel[1]=1 */ + Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres][1]=1(V4=1) Tresult[nres][2]=0(V3=0) */ + Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */ + Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */ + printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4); + k4++;; + } else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */ + k3q= resultmodel[k1]; /* resultmodel[1(V5)] = 25.1=k3q */ + k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ + Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */ + Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */ + Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ + printf("Decoderesult Quantitative nres=%d, V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); + k4q++;; + } + } + + TKresult[nres]=++k; /* Combination for the nresult and the model */ + return (0); +} + +int decodemodel( char model[], int lastobs) + /**< This routine decodes the model and returns: + * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age + * - nagesqr = 1 if age*age in the model, otherwise 0. + * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age + * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age + * - cptcovage number of covariates with age*products =2 + * - cptcovs number of simple covariates + * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10 + * which is a new column after the 9 (ncovcol) variables. + * - if k is a product Vn*Vm, covar[k][i] is filled with correct values for each individual + * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage + * Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6. + * - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 . + */ +/* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */ +{ + int i, j, k, ks, v; + int j1, k1, k2, k3, k4; + char modelsav[80]; + char stra[80], strb[80], strc[80], strd[80],stre[80]; + char *strpt; + + /*removespace(model);*/ + if (strlen(model) >1){ /* If there is at least 1 covariate */ + j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0; + if (strstr(model,"AGE") !=0){ + printf("Error. AGE must be in lower case 'age' model=1+age+%s. ",model); + fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s. ",model);fflush(ficlog); + return 1; + } + if (strstr(model,"v") !=0){ + printf("Error. 'v' must be in upper case 'V' model=%s ",model); + fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog); + return 1; + } + strcpy(modelsav,model); + if ((strpt=strstr(model,"age*age")) !=0){ + printf(" strpt=%s, model=%s\n",strpt, model); + if(strpt != model){ + printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ + corresponding column of parameters.\n",model); + fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ + corresponding column of parameters.\n",model); fflush(ficlog); + return 1; + } + nagesqr=1; + if (strstr(model,"+age*age") !=0) + substrchaine(modelsav, model, "+age*age"); + else if (strstr(model,"age*age+") !=0) + substrchaine(modelsav, model, "age*age+"); + else + substrchaine(modelsav, model, "age*age"); + }else + nagesqr=0; + if (strlen(modelsav) >1){ + j=nbocc(modelsav,'+'); /**< j=Number of '+' */ + j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */ + cptcovs=j+1-j1; /**< Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =5-3=2 */ + cptcovt= j+1; /* Number of total covariates in the model, not including + * cst, age and age*age + * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/ + /* including age products which are counted in cptcovage. + * but the covariates which are products must be treated + * separately: ncovn=4- 2=2 (V1+V3). */ + cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */ + cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ + + + /* Design + * V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight + * < ncovcol=8 > + * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 + * k= 1 2 3 4 5 6 7 8 + * cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8 + * covar[k,i], value of kth covariate if not including age for individual i: + * covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8) + * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[2]=1 Tvar[4]=3 Tvar[8]=8 + * if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and + * Tage[++cptcovage]=k + * if products, new covar are created after ncovcol with k1 + * Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11 + * Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product + * Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8 + * Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2]; + * Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted + * V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 + * < ncovcol=8 > + * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2 + * k= 1 2 3 4 5 6 7 8 9 10 11 12 + * Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8 + * p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} + * p Tprod[1]@2={ 6, 5} + *p Tvard[1][1]@4= {7, 8, 5, 6} + * covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8 + * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; + *How to reorganize? Tvars(orted) + * Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age + * Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} + * {2, 1, 4, 8, 5, 6, 3, 7} + * Struct [] + */ + + /* This loop fills the array Tvar from the string 'model'.*/ + /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */ + /* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */ + /* k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */ + /* k=3 V4 Tvar[k=3]= 4 (from V4) */ + /* k=2 V1 Tvar[k=2]= 1 (from V1) */ + /* k=1 Tvar[1]=2 (from V2) */ + /* k=5 Tvar[5] */ + /* for (k=1; k<=cptcovn;k++) { */ + /* cov[2+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + /* } */ + /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */ + /* + * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ + for(k=cptcovt; k>=1;k--){ /**< Number of covariates not including constant and age, neither age*age*/ + Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0; + } + cptcovage=0; + for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */ + cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right + modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */ /*"V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */ + if (nbocc(modelsav,'+')==0) + strcpy(strb,modelsav); /* and analyzes it */ + /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ + /*scanf("%d",i);*/ + if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V5*age+ V4+V3*age strb=V3*age */ + cutl(strc,strd,strb,'*'); /**< k=1 strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ + if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ + /* covar is not filled and then is empty */ + cptcovprod--; + cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ + Tvar[k]=atoi(stre); /* V2+V1+V5*age+V4+V3*age Tvar[5]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ + Typevar[k]=1; /* 1 for age product */ + cptcovage++; /* Counts the number of covariates which include age as a product */ + Tage[cptcovage]=k; /* V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ + /*printf("stre=%s ", stre);*/ + } else if (strcmp(strd,"age")==0) { /* or age*Vn */ + cptcovprod--; + cutl(stre,strb,strc,'V'); + Tvar[k]=atoi(stre); + Typevar[k]=1; /* 1 for age product */ + cptcovage++; + Tage[cptcovage]=k; + } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ + /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ + cptcovn++; + cptcovprodnoage++;k1++; + cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ + Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but + because this model-covariate is a construction we invent a new column + which is after existing variables ncovcol+nqv+ntv+nqtv + k1 + If already ncovcol=4 and model=V2 + V1 +V1*V4 +age*V3 +V3*V2 + thus after V4 we invent V5 and V6 because age*V3 will be computed in 4 + Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */ + Typevar[k]=2; /* 2 for double fixed dummy covariates */ + cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ + Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ + Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */ + Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ + Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ + k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */ + /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */ + /* Tvar[cptcovt+k2+1]=Tvard[k1][2]; /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */ + /*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */ + /* 1 2 3 4 5 | Tvar[5+1)=1, Tvar[7]=2 */ + for (i=1; i<=lastobs;i++){ + /* Computes the new covariate which is a product of + covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ + covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; + } + } /* End age is not in the model */ + } /* End if model includes a product */ + else { /* not a product */ + /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ + /* scanf("%d",i);*/ + cutl(strd,strc,strb,'V'); + ks++; /**< Number of simple covariates dummy or quantitative, fixe or varying */ + cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */ + Tvar[k]=atoi(strd); + Typevar[k]=0; /* 0 for simple covariates */ + } + strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ + /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); + scanf("%d",i);*/ + } /* end of loop + on total covariates */ + } /* end if strlen(modelsave == 0) age*age might exist */ + } /* end if strlen(model == 0) */ + + /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. + If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ + + /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); + printf("cptcovprod=%d ", cptcovprod); + fprintf(ficlog,"cptcovprod=%d ", cptcovprod); + scanf("%d ",i);*/ + + +/* Until here, decodemodel knows only the grammar (simple, product, age*) of the model but not what kind + of variable (dummy vs quantitative, fixed vs time varying) is behind. But we know the # of each. */ +/* ncovcol= 1, nqv=1 | ntv=2, nqtv= 1 = 5 possible variables data: 2 fixed 3, varying + model= V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place + k = 1 2 3 4 5 6 7 8 9 + Tvar[k]= 5 4 3 1+1+2+1+1=6 5 2 7 1 5 + Typevar[k]= 0 0 0 2 1 0 2 1 0 + Fixed[k] 1 1 1 1 3 0 0 or 2 2 3 + Dummy[k] 1 0 0 0 3 1 1 2 3 + Tmodelind[combination of covar]=k; +*/ +/* Dispatching between quantitative and time varying covariates */ + /* If Tvar[k] >ncovcol it is a product */ + /* Tvar[k] is the value n of Vn with n varying for 1 to nvcol, or p Vp=Vn*Vm for product */ + /* Computing effective variables, ie used by the model, that is from the cptcovt variables */ + printf("Model=1+age+%s\n\ +Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ +Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ +Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); + fprintf(ficlog,"Model=1+age+%s\n\ +Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ +Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ +Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); + for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} + for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ + if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ + Fixed[k]= 0; + Dummy[k]= 0; + ncoveff++; + ncovf++; + nsd++; + modell[k].maintype= FTYPE; + TvarsD[nsd]=Tvar[k]; + TvarsDind[nsd]=k; + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + }else if( Tvar[k] <=ncovcol && Typevar[k]==2){ /* Product of fixed dummy (<=ncovcol) covariates */ + Fixed[k]= 0; + Dummy[k]= 0; + ncoveff++; + ncovf++; + modell[k].maintype= FTYPE; + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */ + Fixed[k]= 0; + Dummy[k]= 1; + nqfveff++; + modell[k].maintype= FTYPE; + modell[k].subtype= FQ; + nsq++; + TvarsQ[nsq]=Tvar[k]; + TvarsQind[nsq]=k; + ncovf++; + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ + TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ + }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */ + Fixed[k]= 1; + Dummy[k]= 0; + ntveff++; /* Only simple time varying dummy variable */ + modell[k].maintype= VTYPE; + modell[k].subtype= VD; + nsd++; + TvarsD[nsd]=Tvar[k]; + TvarsDind[nsd]=k; + ncovv++; /* Only simple time varying variables */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; /* TvarVind[2]=2 TvarVind[3]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ + TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ + TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ + printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv); + printf("Quasi TmodelInvind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv); + }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){ /* Only simple time varying quantitative variable V5*/ + Fixed[k]= 1; + Dummy[k]= 1; + nqtveff++; + modell[k].maintype= VTYPE; + modell[k].subtype= VQ; + ncovv++; /* Only simple time varying variables */ + nsq++; + TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */ + TvarsQind[nsq]=k; + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ + TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ + TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ + TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */ + /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */ + printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv); + printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv); + }else if (Typevar[k] == 1) { /* product with age */ + ncova++; + TvarA[ncova]=Tvar[k]; + TvarAind[ncova]=k; + if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */ + Fixed[k]= 2; + Dummy[k]= 2; + modell[k].maintype= ATYPE; + modell[k].subtype= APFD; + /* ncoveff++; */ + }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/ + Fixed[k]= 2; + Dummy[k]= 3; + modell[k].maintype= ATYPE; + modell[k].subtype= APFQ; /* Product age * fixed quantitative */ + /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */ + }else if( Tvar[k] <=ncovcol+nqv+ntv ){ + Fixed[k]= 3; + Dummy[k]= 2; + modell[k].maintype= ATYPE; + modell[k].subtype= APVD; /* Product age * varying dummy */ + /* ntveff++; /\* Only simple time varying dummy variable *\/ */ + }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 3; + Dummy[k]= 3; + modell[k].maintype= ATYPE; + modell[k].subtype= APVQ; /* Product age * varying quantitative */ + /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */ + } + }else if (Typevar[k] == 2) { /* product without age */ + k1=Tposprod[k]; + if(Tvard[k1][1] <=ncovcol){ + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */ + ncovf++; /* Fixed variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 0; /* or 2 ?*/ + Dummy[k]= 1; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */ + ncovf++; /* Varying variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }else if(Tvard[k1][1] <=ncovcol+nqv){ + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 0; /* or 2 ?*/ + Dummy[k]= 1; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */ + ncovf++; /* Fixed variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }else{ + printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); + fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); + } /*end k1*/ + }else{ + printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); + fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); + } + printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]); + printf(" modell[%d].maintype=%d, modell[%d].subtype=%d\n",k,modell[k].maintype,k,modell[k].subtype); + fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]); + } + /* Searching for doublons in the model */ + for(k1=1; k1<= cptcovt;k1++){ + for(k2=1; k2 Result files
\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):
"); - fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year
",p[1],p[2],agegomp); - for (i=1;i<=2;i++) - fprintf(fichtm," p[%d] = %lf [%f ; %f]
\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); - fprintf(fichtm,"![](\"graphmort.png\")
");
- fprintf(fichtm,"Parameter files
\n\ + - Parameter file: %s.%s
\n\ - Copy of the parameter file: o%s
\n\ - Log file of the run: %s
\n\ - Gnuplot file name: %s
\n\ - Date and time at start: %s- Tolerance for the convergence of the likelihood: ftol=%g \n
- Interval for the elementary matrix (in month): stepm=%d",\ + ftol, stepm); + fprintf(fichtm,"\n
- Number of fixed dummy covariates: ncovcol=%d ", ncovcol); + ncurrv=1; + for(i=ncurrv; i <=ncovcol; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
- Number of fixed quantitative variables: nqv=%d ", nqv); + ncurrv=i; + for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
- Number of time varying (wave varying) dummy covariates: ntv=%d ", ntv); + ncurrv=i; + for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
- Number of time varying quantitative covariates: nqtv=%d ", nqtv); + ncurrv=i; + for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
- Weights column \n
Number of alive states: nlstate=%d
Number of death states (not really implemented): ndeath=%d \n - Number of waves: maxwav=%d \n
- Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n
- Does the weight column be taken into account (1), or not (0): weight=%d
- model=1+age+%s\n \ +
"); + + jj1=0; + + fprintf(fichtm," \n
- ");
+ for(nres=1; nres <= nresult; nres++) /* For each resultline */
+ for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
+ if(m != 1 && TKresult[nres]!= k1)
+ continue;
+ jj1++;
+ if (cptcovn > 0) {
+ fprintf(fichtm,"\n
"); + jj1=0; - m=cptcoveff; - if (cptcovn < 1) {m=1;ncodemax[1]=1;} + for(nres=1; nres <= nresult; nres++) /* For each resultline */ + for(k1=1; k1<=m;k1++){ /* For each combination of covariate */ + if(m != 1 && TKresult[nres]!= k1) + continue; - jj1=0; - for(k1=1; k1<=m;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ + /* for(i1=1; i1<=ncodemax[k1];i1++){ */ jj1++; if (cptcovn > 0) { + fprintf(fichtm,"\n
");
+
fprintf(fichtm,"
************ Results for covariates");
- for (cpt=1; cpt<=cptcoveff;cpt++)
- fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
- fprintf(fichtm," ************\n
");
+ for (cpt=1; cpt<=cptcoveff;cpt++){
+ fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
+ printf(" V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);fflush(stdout);
+ /* 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); */
+ }
+ for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+ fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);fflush(stdout);
+ }
+
+ /* if(nqfveff+nqtveff 0) */ /* Test to be done */
+ fprintf(fichtm," (model=%s) ************\n
",model);
+ if(invalidvarcomb[k1]){
+ fprintf(fichtm,"\nCombination (%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: logit(pij)=log(pij/pii)= aij+ bij age+%s) as a function of age: %s_%d-1-%d.svg
\
+![](\"%s_%d-1-%d.svg\")
",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Pij */
- fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png
\
-![](\"%s%d1.png\")
",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ fprintf(fichtm,"
\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s_%d-2-%d.svg
\
+![](\"%s_%d-2-%d.svg\")
",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Quasi-incidences */
- fprintf(fichtm,"
- Pij 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: %s%d2.png
\
-![](\"%s%d2.png\")
",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
- /* Stable prevalence in each health state */
- for(cpt=1; cpt
\
-![](\"%s%d%d.png\")
",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+ 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-%d.svg
\
+![](\"%s_%d-3-%d.svg\")
",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
+ /* Survival functions (period) in state j */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\
+![](\"%s_%d-%d-%d.svg\")
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+ }
+ /* State specific survival functions (period) */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Survival functions in state %d and in any other live state (total).\
+ And probability to be observed in various states (up to %d) being in state %d at different ages. \
+ %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+ }
+ /* Period (forward 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 for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
\
+", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+ }
+ if(prevbcast==1){
+ /* Backward prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. %s_%d-%d-%d.svg
\
+", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
+ }
+ }
+ if(prevfcast==1){
+ /* Projection of prevalence up to period (forward 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 and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg", dateprev1, dateprev2, mobilavproj, dateprojd, dateprojf, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"F_"),subdirf2(optionfilefiname,"F_"));
+ fprintf(fichtm,"",
+ subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
}
+ }
+ if(prevbcast==1){
+ /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
+ from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \
+ account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
+with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));
+ fprintf(fichtm," ", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ }
+ }
+
for(cpt=1; cpt<=nlstate;cpt++) {
- fprintf(fichtm,"\n
- Health life expectancies by age and initial health state (%d): %s%d%d.png
\
-",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+ fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d-%d-%d.svg",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"E_"),subdirf2(optionfilefiname,"E_"));
+ fprintf(fichtm,"", subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres );
}
- } /* end i1 */
- }/* End k1 */
- fprintf(fichtm,"
\n", + estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); + fprintf(fichtm,"\ + - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the forward (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,"\ + - 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,"\ + - Standard deviation of forward (period) prevalences: %s
\n",\ + subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); /* if(popforecast==1) fprintf(fichtm,"\n */ /* - Prevalences forecasting: f%s
\n */ @@ -3233,42 +7430,118 @@ fprintf(fichtm," \n
"); + fflush(fichtm); + + m=pow(2,cptcoveff); + if (cptcovn < 1) {m=1;ncodemax[1]=1;} + + fprintf(fichtm,"
"); + + jj1=0; + + fprintf(fichtm," \n
- ");
+ for(nres=1; nres <= nresult; nres++) /* For each resultline */
+ for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
+ if(m != 1 && TKresult[nres]!= k1)
+ continue;
+ jj1++;
+ if (cptcovn > 0) {
+ fprintf(fichtm,"\n
");
+
fprintf(fichtm,"
************ Results for covariates");
- for (cpt=1; cpt<=cptcoveff;cpt++)
- fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
- fprintf(fichtm," ************\n
");
+ for (cpt=1; cpt<=cptcoveff;cpt++){ /**< cptcoveff number of variables */
+ fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);
+ printf(" V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);fflush(stdout);
+ /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+ }
+ for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+ fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ }
+
+ fprintf(fichtm," (model=%s) ************\n
",model);
+
+ if(invalidvarcomb[k1]){
+ fprintf(fichtm,"\nCombination (%d) ignored because no cases
\n",k1);
+ continue;
+ }
}
for(cpt=1; cpt<=nlstate;cpt++) {
- fprintf(fichtm,"
- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): %s%d%d.png
\
-",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
+ fprintf(fichtm,"\n
- Observed (cross-sectional with mov_average=%d) and period (incidence based) \
+prevalence (with 95%% confidence interval) in state (%d): %s_%d-%d-%d.svg",mobilav,cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s)\n
",subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
+ fprintf(fichtm,"",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
}
fprintf(fichtm,"\n
- Total life expectancy by age and \
-health expectancies in states (1) and (2): %s%d.png
\
-![](\"%s%d.png\")
",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
- } /* end i1 */
- }/* End k1 */
- fprintf(fichtm,"
",subdirf2(optionfilefiname,"T_"),subdirf2(optionfilefiname,"T_")); + fprintf(fichtm,"
![](\"%s_%d-%d.svg\")
",subdirf2(optionfilefiname,"E_"),k1,nres);
+ /* } /\* end i1 *\/ */
+ }/* End k1 */
+ }/* End nres */
+ fprintf(fichtm,"%s \ + fprintf(fichtmcov,"\n
%s \
\n\ -Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n",\ - fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ + optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(fichtm,"\n
%s \ + fprintf(fichtm,"\n\n\n
\nSponsored by Copyright (C) 2002-2015 INED-EUROREVES-Institut de longeÌviteÌ-2013-2016-Japan Society for the Promotion of Sciences 日本å¦è¡“振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015-2018
\ +
\n\ +IMaCh-%s
%s \
\n\ -Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\ \n\
\
Total number of observations=%d
\n\ + fprintf(fichtm,"
Parameter line 2
Diagram of states %s_.svg
\n\ +![](\"%s_.svg\")
", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"));
+
+
+ fprintf(fichtm,"\nSome descriptive statistics
\nNumber of (used) observations=%d
\n\ Youngest age at first (selected) pass %.2f, oldest age %.2f
\n\ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ - imx,agemin,agemax,jmin,jmax,jmean); + imx,agemin,agemax,jmin,jmax,jmean); pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ - oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ - newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ - savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ - oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ - - + oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ + newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ + savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ + oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ + /* For Powell, parameters are in a vector p[] starting at p[1] so we point p on param[1][1] so that p[1] maps on param[1][1][1] */ p=param[1][1]; /* *(*(*(param +1)+1)+0) */ globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ + /* For mortality only */ if (mle==-3){ - ximort=matrix(1,NDIM,1,NDIM); - cens=ivector(1,n); - ageexmed=vector(1,n); - agecens=vector(1,n); - dcwave=ivector(1,n); - + ximort=matrix(1,NDIM,1,NDIM); + for(i=1;i<=NDIM;i++) + for(j=1;j<=NDIM;j++) + ximort[i][j]=0.; + /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ + cens=ivector(firstobs,lastobs); + ageexmed=vector(firstobs,lastobs); + agecens=vector(firstobs,lastobs); + dcwave=ivector(firstobs,lastobs); + for (i=1; i<=imx; i++){ dcwave[i]=-1; - for (j=1; j<=lastpass; j++) - if (s[j][i]>nlstate) { - dcwave[i]=j; + for (m=firstpass; m<=lastpass; m++) + if (s[m][i]>nlstate) { + dcwave[i]=m; /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ break; } } - + for (i=1; i<=imx; i++) { if (wav[i]>0){ ageexmed[i]=agev[mw[1][i]][i]; - j=wav[i];agecens[i]=1.; - if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i]; - cens[i]=1; + j=wav[i]; + agecens[i]=1.; - if (ageexmed[i]<1) cens[i]=-1; - if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ; + if (ageexmed[i]> 1 && wav[i] > 0){ + agecens[i]=agev[mw[j][i]][i]; + cens[i]= 1; + }else if (ageexmed[i]< 1) + cens[i]= -1; + if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass) + cens[i]=0 ; } else cens[i]=-1; } @@ -4795,100 +12329,354 @@ Interval (in months) between two waves: for (j=1;j<=NDIM;j++) ximort[i][j]=(i == j ? 1.0 : 0.0); } + + p[1]=0.0268; p[NDIM]=0.083; + /* printf("%lf %lf", p[1], p[2]); */ + + +#ifdef GSL + printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); +#else + printf("Powell\n"); fprintf(ficlog,"Powell\n"); +#endif + strcpy(filerespow,"POW-MORT_"); + strcat(filerespow,fileresu); + if((ficrespow=fopen(filerespow,"w"))==NULL) { + printf("Problem with resultfile: %s\n", filerespow); + fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); + } +#ifdef GSL + fprintf(ficrespow,"# GSL optimization\n# iter -2*LL"); +#else + fprintf(ficrespow,"# Powell\n# iter -2*LL"); +#endif + /* for (i=1;i<=nlstate;i++) + for(j=1;j<=nlstate+ndeath;j++) + if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); + */ + fprintf(ficrespow,"\n"); +#ifdef GSL + /* gsl starts here */ + T = gsl_multimin_fminimizer_nmsimplex; + gsl_multimin_fminimizer *sfm = NULL; + gsl_vector *ss, *x; + gsl_multimin_function minex_func; - p[1]=0.1; p[2]=0.1; - /*printf("%lf %lf", p[1], p[2]);*/ + /* Initial vertex size vector */ + ss = gsl_vector_alloc (NDIM); + if (ss == NULL){ + GSL_ERROR_VAL ("failed to allocate space for ss", GSL_ENOMEM, 0); + } + /* Set all step sizes to 1 */ + gsl_vector_set_all (ss, 0.001); + + /* Starting point */ - printf("Powell\n"); fprintf(ficlog,"Powell\n"); - strcpy(filerespow,"pow-mort"); - strcat(filerespow,fileres); - if((ficrespow=fopen(filerespow,"w"))==NULL) { - printf("Problem with resultfile: %s\n", filerespow); - fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); - } - fprintf(ficrespow,"# Powell\n# iter -2*LL"); - /* for (i=1;i<=nlstate;i++) - for(j=1;j<=nlstate+ndeath;j++) - if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); - */ - fprintf(ficrespow,"\n"); + x = gsl_vector_alloc (NDIM); + + if (x == NULL){ + gsl_vector_free(ss); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + /* Initialize method and iterate */ + /* p[1]=0.0268; p[NDIM]=0.083; */ + /* gsl_vector_set(x, 0, 0.0268); */ + /* gsl_vector_set(x, 1, 0.083); */ + gsl_vector_set(x, 0, p[1]); + gsl_vector_set(x, 1, p[2]); + + minex_func.f = &gompertz_f; + minex_func.n = NDIM; + minex_func.params = (void *)&p; /* ??? */ + + sfm = gsl_multimin_fminimizer_alloc (T, NDIM); + gsl_multimin_fminimizer_set (sfm, &minex_func, x, ss); + + printf("Iterations beginning .....\n\n"); + printf("Iter. # Intercept Slope -Log Likelihood Simplex size\n"); - powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); + iteri=0; + while (rval == GSL_CONTINUE){ + iteri++; + status = gsl_multimin_fminimizer_iterate(sfm); + + if (status) printf("error: %s\n", gsl_strerror (status)); + fflush(0); + + if (status) + break; + + rval = gsl_multimin_test_size (gsl_multimin_fminimizer_size (sfm), 1e-6); + ssval = gsl_multimin_fminimizer_size (sfm); + + if (rval == GSL_SUCCESS) + printf ("converged to a local maximum at\n"); + + printf("%5d ", iteri); + for (it = 0; it < NDIM; it++){ + printf ("%10.5f ", gsl_vector_get (sfm->x, it)); + } + printf("f() = %-10.5f ssize = %.7f\n", sfm->fval, ssval); + } + + printf("\n\n Please note: Program should be run many times with varying starting points to detemine global maximum\n\n"); + + gsl_vector_free(x); /* initial values */ + gsl_vector_free(ss); /* inital step size */ + for (it=0; it
| Model= | 1 | + age | "); + if(nagesqr==1){ + printf(" + age*age "); + fprintf(ficres," + age*age "); + fprintf(ficlog," + age*age "); + fprintf(fichtm, "+ age*age | "); + } + for(j=1;j <=ncovmodel-2;j++){ + if(Typevar[j]==0) { + printf(" + V%d ",Tvar[j]); + fprintf(ficres," + V%d ",Tvar[j]); + fprintf(ficlog," + V%d ",Tvar[j]); + fprintf(fichtm, "+ V%d | ",Tvar[j]); + }else if(Typevar[j]==1) { + printf(" + V%d*age ",Tvar[j]); + fprintf(ficres," + V%d*age ",Tvar[j]); + fprintf(ficlog," + V%d*age ",Tvar[j]); + fprintf(fichtm, "+ V%d*age | ",Tvar[j]); + }else if(Typevar[j]==2) { + printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficres," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(fichtm, "+ V%d*V%d | ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + } + } + printf("\n"); + fprintf(ficres,"\n"); + fprintf(ficlog,"\n"); + fprintf(fichtm, "
|---|---|---|---|---|---|---|
| %1d%1d | ",i,k); for(j=1; j <=ncovmodel; j++){ - printf("%f ",p[jk]); - fprintf(ficlog,"%f ",p[jk]); - fprintf(ficres,"%f ",p[jk]); + printf("%12.7f ",p[jk]); + fprintf(ficlog,"%12.7f ",p[jk]); + fprintf(ficres,"%12.7f ",p[jk]); + fprintf(fichtm, "%12.7f | ",p[jk]); jk++; } printf("\n"); fprintf(ficlog,"\n"); fprintf(ficres,"\n"); + fprintf(fichtm, "
The Wald test results are output only if the maximimzation of the Likelihood is performed (mle=1)\nParameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities and its graphs'.\n",optionfilehtmcov); + fprintf(fichtm,"\n
| Model= | 1 | + age | "); + if(nagesqr==1){ + printf(" + age*age "); + fprintf(ficres," + age*age "); + fprintf(ficlog," + age*age "); + fprintf(fichtm, "+ age*age | "); + } + for(j=1;j <=ncovmodel-2;j++){ + if(Typevar[j]==0) { + printf(" + V%d ",Tvar[j]); + fprintf(fichtm, "+ V%d | ",Tvar[j]); + }else if(Typevar[j]==1) { + printf(" + V%d*age ",Tvar[j]); + fprintf(fichtm, "+ V%d*age | ",Tvar[j]); + }else if(Typevar[j]==2) { + fprintf(fichtm, "+ V%d*V%d | ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + } + } + fprintf(fichtm, "
|---|---|---|---|---|---|---|
| %1d%1d | ",i,k); + for(j=1; j <=ncovmodel; j++){ + wald=p[jk]/sqrt(matcov[jk][jk]); + printf("%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); + fprintf(ficlog,"%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); + if(fabs(wald) > 1.96){ + fprintf(fichtm, "%12.7f (%12.7f)",p[jk],sqrt(matcov[jk][jk])); + }else{ + fprintf(fichtm, " | %12.7f (%12.7f)",p[jk],sqrt(matcov[jk][jk])); + } + fprintf(fichtm,"W=%8.3f",wald); + fprintf(fichtm,"[%12.7f;%12.7f] | ", p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); + jk++; + } + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(fichtm, "
Local time at start %s
Local time at end %s
",strstart, strtend); + fprintf(fichtm,"
Local time at start %s
Local time at end %s
\n",strstart, strtend); fclose(fichtm); + fprintf(fichtmcov,"
Local time at start %s
Local time at end %s
\n",strstart, strtend); fclose(fichtmcov); fclose(ficgp); fclose(ficlog); /*------ End -----------*/ + - chdir(path); - strcpy(plotcmd,"\""); - strcat(plotcmd,pathimach); - strcat(plotcmd,GNUPLOTPROGRAM); - strcat(plotcmd,"\""); - strcat(plotcmd," "); - strcat(plotcmd,optionfilegnuplot); - printf("Starting graphs with: %s",plotcmd);fflush(stdout); +/* Executes gnuplot */ + + printf("Before Current directory %s!\n",pathcd); +#ifdef WIN32 + if (_chdir(pathcd) != 0) + printf("Can't move to directory %s!\n",path); + if(_getcwd(pathcd,MAXLINE) > 0) +#else + if(chdir(pathcd) != 0) + printf("Can't move to directory %s!\n", path); + if (getcwd(pathcd, MAXLINE) > 0) +#endif + printf("Current directory %s!\n",pathcd); + /*strcat(plotcmd,CHARSEPARATOR);*/ + sprintf(plotcmd,"gnuplot"); +#ifdef _WIN32 + sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); +#endif + if(!stat(plotcmd,&info)){ + printf("Error or gnuplot program not found: '%s'\n",plotcmd);fflush(stdout); + if(!stat(getenv("GNUPLOTBIN"),&info)){ + printf("Error or gnuplot program not found: '%s' Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout); + }else + strcpy(pplotcmd,plotcmd); +#ifdef __unix + strcpy(plotcmd,GNUPLOTPROGRAM); + if(!stat(plotcmd,&info)){ + printf("Error gnuplot program not found: '%s'\n",plotcmd);fflush(stdout); + }else + strcpy(pplotcmd,plotcmd); +#endif + }else + strcpy(pplotcmd,plotcmd); + + sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); + printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); + strcpy(pplotcmd,plotcmd); + if((outcmd=system(plotcmd)) != 0){ - printf(" Problem with gnuplot\n"); + printf("Error in gnuplot, command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); + printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); + sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); + if((outcmd=system(plotcmd)) != 0){ + printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); + strcpy(plotcmd,pplotcmd); + } } - printf(" Wait..."); + printf(" Successful, please wait..."); while (z[0] != 'q') { /* chdir(path); */ - printf("\nType e to edit output files, g to graph again and q for exiting: "); + printf("\nType e to edit results with your browser, g to graph again and q for exit: "); scanf("%s",z); /* if (z[0] == 'c') system("./imach"); */ if (z[0] == 'e') { - printf("Starting browser with: %s",optionfilehtm);fflush(stdout); - system(optionfilehtm); +#ifdef __APPLE__ + sprintf(pplotcmd, "open %s", optionfilehtm); +#elif __linux + sprintf(pplotcmd, "xdg-open %s", optionfilehtm); +#else + sprintf(pplotcmd, "%s", optionfilehtm); +#endif + printf("Starting browser with: %s",pplotcmd);fflush(stdout); + system(pplotcmd); } else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } - end: +end: while (z[0] != 'q') { - printf("\nType q for exiting: "); + printf("\nType q for exiting: "); fflush(stdout); scanf("%s",z); } + printf("End\n"); + exit(0); } - - -