|
version 1.34, 2002/03/13 17:19:16
|
version 1.44, 2002/05/24 13:01:48
|
|
Line 14
|
Line 14
|
| model. More health states you consider, more time is necessary to reach the
|
model. More health states you consider, more time is necessary to reach the
|
| Maximum Likelihood of the parameters involved in the model. The
|
Maximum Likelihood of the parameters involved in the model. The
|
| simplest model is the multinomial logistic model where pij is the
|
simplest model is the multinomial logistic model where pij is the
|
| probabibility to be observed in state j at the second wave
|
probability to be observed in state j at the second wave
|
| conditional to be observed in state i at the first wave. Therefore
|
conditional to be observed in state i at the first wave. Therefore
|
| the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
|
the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
|
| 'age' is age and 'sex' is a covariate. If you want to have a more
|
'age' is age and 'sex' is a covariate. If you want to have a more
|
|
Line 56
|
Line 56
|
| #include <unistd.h>
|
#include <unistd.h>
|
|
|
|
| #define MAXLINE 256
|
#define MAXLINE 256
|
| #define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"
|
#define GNUPLOTPROGRAM "gnuplot"
|
| |
/*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
|
| #define FILENAMELENGTH 80
|
#define FILENAMELENGTH 80
|
| /*#define DEBUG*/
|
/*#define DEBUG*/
|
| #define windows
|
#define windows
|
|
Line 135 int imx;
|
Line 136 int imx;
|
| int stepm;
|
int stepm;
|
| /* Stepm, step in month: minimum step interpolation*/
|
/* Stepm, step in month: minimum step interpolation*/
|
|
|
|
| |
int estepm;
|
| |
/* Estepm, step in month to interpolate survival function in order to approximate Life Expectancy*/
|
| |
|
| int m,nb;
|
int m,nb;
|
| int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
|
int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
|
| double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
|
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
|
|
Line 686 double **prevalim(double **prlim, int nl
|
Line 690 double **prevalim(double **prlim, int nl
|
|
|
|
| for (k=1; k<=cptcovn;k++) {
|
for (k=1; k<=cptcovn;k++) {
|
| cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
|
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
|
| /*printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);*/
|
/* printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/
|
| }
|
}
|
| for (k=1; k<=cptcovage;k++)
|
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
|
| cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
|
|
| for (k=1; k<=cptcovprod;k++)
|
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]]];
|
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
|
|
|
|
| /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
|
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
|
| /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
|
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
|
|
|
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
|
| out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
|
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
|
|
|
|
| savm=oldm;
|
savm=oldm;
|
|
Line 1178 void lubksb(double **a, int n, int *indx
|
Line 1181 void lubksb(double **a, int n, int *indx
|
| /************ Frequencies ********************/
|
/************ Frequencies ********************/
|
| void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
|
void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
|
| { /* Some frequencies */
|
{ /* Some frequencies */
|
|
|
|
| int i, m, jk, k1,i1, j1, bool, z1,z2,j;
|
int i, m, jk, k1,i1, j1, bool, z1,z2,j;
|
| double ***freq; /* Frequencies */
|
double ***freq; /* Frequencies */
|
| double *pp;
|
double *pp;
|
| double pos, k2, dateintsum=0,k2cpt=0;
|
double pos, k2, dateintsum=0,k2cpt=0;
|
| FILE *ficresp;
|
FILE *ficresp;
|
| char fileresp[FILENAMELENGTH];
|
char fileresp[FILENAMELENGTH];
|
|
|
|
| pp=vector(1,nlstate);
|
pp=vector(1,nlstate);
|
| probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
| strcpy(fileresp,"p");
|
strcpy(fileresp,"p");
|
|
Line 1196 void freqsummary(char fileres[], int ag
|
Line 1199 void freqsummary(char fileres[], int ag
|
| }
|
}
|
| freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
|
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
|
| j1=0;
|
j1=0;
|
|
|
|
| j=cptcoveff;
|
j=cptcoveff;
|
| if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
|
|
|
| for(k1=1; k1<=j;k1++){
|
for(k1=1; k1<=j;k1++){
|
| for(i1=1; i1<=ncodemax[k1];i1++){
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
| j1++;
|
j1++;
|
| /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
|
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
|
| scanf("%d", i);*/
|
scanf("%d", i);*/
|
| for (i=-1; i<=nlstate+ndeath; i++)
|
for (i=-1; i<=nlstate+ndeath; i++)
|
| for (jk=-1; jk<=nlstate+ndeath; jk++)
|
for (jk=-1; jk<=nlstate+ndeath; jk++)
|
| for(m=agemin; m <= agemax+3; m++)
|
for(m=agemin; m <= agemax+3; m++)
|
| freq[i][jk][m]=0;
|
freq[i][jk][m]=0;
|
|
|
|
| dateintsum=0;
|
dateintsum=0;
|
| k2cpt=0;
|
k2cpt=0;
|
| for (i=1; i<=imx; i++) {
|
for (i=1; i<=imx; i++) {
|
| bool=1;
|
bool=1;
|
| if (cptcovn>0) {
|
if (cptcovn>0) {
|
| for (z1=1; z1<=cptcoveff; z1++)
|
for (z1=1; z1<=cptcoveff; z1++)
|
| if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
|
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
|
| bool=0;
|
bool=0;
|
| }
|
}
|
| if (bool==1) {
|
if (bool==1) {
|
| for(m=firstpass; m<=lastpass; m++){
|
for(m=firstpass; m<=lastpass; m++){
|
| k2=anint[m][i]+(mint[m][i]/12.);
|
k2=anint[m][i]+(mint[m][i]/12.);
|
| if ((k2>=dateprev1) && (k2<=dateprev2)) {
|
if ((k2>=dateprev1) && (k2<=dateprev2)) {
|
| if(agev[m][i]==0) agev[m][i]=agemax+1;
|
if(agev[m][i]==0) agev[m][i]=agemax+1;
|
| if(agev[m][i]==1) agev[m][i]=agemax+2;
|
if(agev[m][i]==1) agev[m][i]=agemax+2;
|
| freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
|
if (m<lastpass) {
|
| freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
|
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
|
| if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
|
freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
|
| dateintsum=dateintsum+k2;
|
}
|
| k2cpt++;
|
|
| }
|
if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
|
|
|
dateintsum=dateintsum+k2;
|
| }
|
k2cpt++;
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
| |
}
|
| |
}
|
|
|
|
| fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
|
|
|
| if (cptcovn>0) {
|
if (cptcovn>0) {
|
| fprintf(ficresp, "\n#********** Variable ");
|
fprintf(ficresp, "\n#********** Variable ");
|
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
| fprintf(ficresp, "**********\n#");
|
fprintf(ficresp, "**********\n#");
|
| }
|
}
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
|
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
|
| fprintf(ficresp, "\n");
|
fprintf(ficresp, "\n");
|
|
|
|
| for(i=(int)agemin; i <= (int)agemax+3; i++){
|
for(i=(int)agemin; i <= (int)agemax+3; i++){
|
| if(i==(int)agemax+3)
|
if(i==(int)agemax+3)
|
| printf("Total");
|
printf("Total");
|
| else
|
else
|
| printf("Age %d", i);
|
printf("Age %d", i);
|
| for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
| for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
|
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
|
| pp[jk] += freq[jk][m][i];
|
pp[jk] += freq[jk][m][i];
|
| }
|
}
|
| for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
| for(m=-1, pos=0; m <=0 ; m++)
|
for(m=-1, pos=0; m <=0 ; m++)
|
| pos += freq[jk][m][i];
|
pos += freq[jk][m][i];
|
| if(pp[jk]>=1.e-10)
|
if(pp[jk]>=1.e-10)
|
| printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
|
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
|
| else
|
else
|
| printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
|
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
|
| }
|
}
|
|
|
|
| for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
| for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
|
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
|
| pp[jk] += freq[jk][m][i];
|
pp[jk] += freq[jk][m][i];
|
| }
|
}
|
|
|
|
| for(jk=1,pos=0; jk <=nlstate ; jk++)
|
for(jk=1,pos=0; jk <=nlstate ; jk++)
|
| pos += pp[jk];
|
pos += pp[jk];
|
| for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
| if(pos>=1.e-5)
|
if(pos>=1.e-5)
|
| printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
|
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
|
| else
|
else
|
| printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
|
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
|
| if( i <= (int) agemax){
|
if( i <= (int) agemax){
|
| if(pos>=1.e-5){
|
if(pos>=1.e-5){
|
| fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
|
fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
|
| probs[i][jk][j1]= pp[jk]/pos;
|
probs[i][jk][j1]= pp[jk]/pos;
|
| /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
|
/*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
|
| |
}
|
| |
else
|
| |
fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
|
| |
}
|
| }
|
}
|
| else
|
|
| fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
|
for(jk=-1; jk <=nlstate+ndeath; jk++)
|
| |
for(m=-1; m <=nlstate+ndeath; m++)
|
| |
if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
|
| |
if(i <= (int) agemax)
|
| |
fprintf(ficresp,"\n");
|
| |
printf("\n");
|
| }
|
}
|
| }
|
}
|
| for(jk=-1; jk <=nlstate+ndeath; jk++)
|
}
|
| for(m=-1; m <=nlstate+ndeath; m++)
|
|
| if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
|
|
| if(i <= (int) agemax)
|
|
| fprintf(ficresp,"\n");
|
|
| printf("\n");
|
|
| }
|
|
| }
|
|
| }
|
|
| dateintmean=dateintsum/k2cpt;
|
dateintmean=dateintsum/k2cpt;
|
|
|
|
| fclose(ficresp);
|
fclose(ficresp);
|
| free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
| free_vector(pp,1,nlstate);
|
free_vector(pp,1,nlstate);
|
|
|
|
| /* End of Freq */
|
/* End of Freq */
|
| }
|
}
|
|
|
|
|
Line 1324 void prevalence(int agemin, float agemax
|
Line 1330 void prevalence(int agemin, float agemax
|
| j=cptcoveff;
|
j=cptcoveff;
|
| if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
|
|
|
| for(k1=1; k1<=j;k1++){
|
for(k1=1; k1<=j;k1++){
|
| for(i1=1; i1<=ncodemax[k1];i1++){
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
| j1++;
|
j1++;
|
|
|
|
| for (i=-1; i<=nlstate+ndeath; i++)
|
for (i=-1; i<=nlstate+ndeath; i++)
|
| for (jk=-1; jk<=nlstate+ndeath; jk++)
|
for (jk=-1; jk<=nlstate+ndeath; jk++)
|
| for(m=agemin; m <= agemax+3; m++)
|
for(m=agemin; m <= agemax+3; m++)
|
|
Line 1346 void prevalence(int agemin, float agemax
|
Line 1352 void prevalence(int agemin, float agemax
|
| if ((k2>=dateprev1) && (k2<=dateprev2)) {
|
if ((k2>=dateprev1) && (k2<=dateprev2)) {
|
| if(agev[m][i]==0) agev[m][i]=agemax+1;
|
if(agev[m][i]==0) agev[m][i]=agemax+1;
|
| if(agev[m][i]==1) agev[m][i]=agemax+2;
|
if(agev[m][i]==1) agev[m][i]=agemax+2;
|
| freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
|
if (m<lastpass) {
|
| /* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i]; */
|
if (calagedate>0)
|
| |
freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
|
| |
else
|
| |
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
|
| |
freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i];
|
| |
}
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
| for(i=(int)agemin; i <= (int)agemax+3; i++){
|
for(i=(int)agemin; i <= (int)agemax+3; i++){
|
| for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
| for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
|
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
|
| pp[jk] += freq[jk][m][i];
|
pp[jk] += freq[jk][m][i];
|
| }
|
}
|
| for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
| for(m=-1, pos=0; m <=0 ; m++)
|
for(m=-1, pos=0; m <=0 ; m++)
|
| pos += freq[jk][m][i];
|
pos += freq[jk][m][i];
|
| }
|
}
|
|
|
|
| for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
| for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
|
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
|
| pp[jk] += freq[jk][m][i];
|
pp[jk] += freq[jk][m][i];
|
| }
|
|
|
|
|
| for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
|
|
|
|
|
| for(jk=1; jk <=nlstate ; jk++){
|
|
| if( i <= (int) agemax){
|
|
| if(pos>=1.e-5){
|
|
| probs[i][jk][j1]= pp[jk]/pos;
|
|
| }
|
|
| }
|
|
| }
|
|
|
|
|
| }
|
}
|
| |
|
| |
for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
|
| |
|
| |
for(jk=1; jk <=nlstate ; jk++){
|
| |
if( i <= (int) agemax){
|
| |
if(pos>=1.e-5){
|
| |
probs[i][jk][j1]= pp[jk]/pos;
|
| |
}
|
| |
}
|
| |
}
|
| |
|
| |
}
|
| }
|
}
|
| }
|
}
|
|
|
|
|
|
|
| free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
| free_vector(pp,1,nlstate);
|
free_vector(pp,1,nlstate);
|
|
Line 1497 void tricode(int *Tvar, int **nbcode, in
|
Line 1508 void tricode(int *Tvar, int **nbcode, in
|
| for (k=0; k<=19; k++) {
|
for (k=0; k<=19; k++) {
|
| if (Ndum[k] != 0) {
|
if (Ndum[k] != 0) {
|
| nbcode[Tvar[j]][ij]=k;
|
nbcode[Tvar[j]][ij]=k;
|
| |
|
| ij++;
|
ij++;
|
| }
|
}
|
| if (ij > ncodemax[j]) break;
|
if (ij > ncodemax[j]) break;
|
|
Line 1524 void tricode(int *Tvar, int **nbcode, in
|
Line 1536 void tricode(int *Tvar, int **nbcode, in
|
|
|
|
| /*********** Health Expectancies ****************/
|
/*********** Health Expectancies ****************/
|
|
|
|
| void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij)
|
void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov )
|
| |
|
| {
|
{
|
| /* Health expectancies */
|
/* Health expectancies */
|
| int i, j, nhstepm, hstepm, h, nstepm, k;
|
int i, j, nhstepm, hstepm, h, nstepm, k, cptj;
|
| double age, agelim,hf;
|
double age, agelim, hf;
|
| double ***p3mat;
|
double ***p3mat,***varhe;
|
| |
double **dnewm,**doldm;
|
| |
double *xp;
|
| |
double **gp, **gm;
|
| |
double ***gradg, ***trgradg;
|
| |
int theta;
|
| |
|
| |
varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage);
|
| |
xp=vector(1,npar);
|
| |
dnewm=matrix(1,nlstate*2,1,npar);
|
| |
doldm=matrix(1,nlstate*2,1,nlstate*2);
|
|
|
|
| fprintf(ficreseij,"# Health expectancies\n");
|
fprintf(ficreseij,"# Health expectancies\n");
|
| fprintf(ficreseij,"# Age");
|
fprintf(ficreseij,"# Age");
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| for(j=1; j<=nlstate;j++)
|
for(j=1; j<=nlstate;j++)
|
| fprintf(ficreseij," %1d-%1d",i,j);
|
fprintf(ficreseij," %1d-%1d (SE)",i,j);
|
| fprintf(ficreseij,"\n");
|
fprintf(ficreseij,"\n");
|
|
|
|
| k=1; /* For example stepm=6 months */
|
if(estepm < stepm){
|
| hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */
|
printf ("Problem %d lower than %d\n",estepm, stepm);
|
| hstepm=stepm; /* or (b) We decided to compute the life expectancy with the smallest unit */
|
}
|
| |
else hstepm=estepm;
|
| |
/* We compute the life expectancy from trapezoids spaced every estepm months
|
| |
* This is mainly to measure the difference between two models: for example
|
| |
* if stepm=24 months pijx are given only every 2 years and by summing them
|
| |
* we are calculating an estimate of the Life Expectancy assuming a linear
|
| |
* progression inbetween and thus overestimating or underestimating according
|
| |
* to the curvature of the survival function. If, for the same date, we
|
| |
* estimate the model with stepm=1 month, we can keep estepm to 24 months
|
| |
* to compare the new estimate of Life expectancy with the same linear
|
| |
* hypothesis. A more precise result, taking into account a more precise
|
| |
* curvature will be obtained if estepm is as small as stepm. */
|
| |
|
| |
/* 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.
|
/* 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
|
nhstepm is the number of hstepm from age to agelim
|
| nstepm is the number of stepm from age to agelin.
|
nstepm is the number of stepm from age to agelin.
|
| Look at hpijx to understand the reason of that which relies in memory size
|
Look at hpijx to understand the reason of that which relies in memory size
|
| and note for a fixed period like k years */
|
and note for a fixed period like estepm months */
|
| /* We decided (b) to get a life expectancy respecting the most precise curvature of the
|
/* 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
|
survival function given by stepm (the optimization length). Unfortunately it
|
| means that if the survival funtion is printed only each two years of age and if
|
means that if the survival funtion is printed only each two years of age and if
|
| you sum them up and add 1 year (area under the trapezoids) you won't get the same
|
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.
|
results. So we changed our mind and took the option of the best precision.
|
| */
|
*/
|
| hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */
|
hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
|
|
|
|
| agelim=AGESUP;
|
agelim=AGESUP;
|
| for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
|
for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
|
|
Line 1562 void evsij(char fileres[], double ***eij
|
Line 1598 void evsij(char fileres[], double ***eij
|
| /* if (stepm >= YEARM) hstepm=1;*/
|
/* if (stepm >= YEARM) hstepm=1;*/
|
| nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
|
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
|
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
| |
gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);
|
| |
gp=matrix(0,nhstepm,1,nlstate*2);
|
| |
gm=matrix(0,nhstepm,1,nlstate*2);
|
| |
|
| /* Computed by stepm unit matrices, product of hstepm matrices, stored
|
/* Computed by stepm unit matrices, product of hstepm matrices, stored
|
| in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
|
in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
|
| hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);
|
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);
|
| |
|
| |
|
| hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
|
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
|
| |
|
| |
/* Computing Variances of health expectancies */
|
| |
|
| |
for(theta=1; theta <=npar; theta++){
|
| |
for(i=1; i<=npar; i++){
|
| |
xp[i] = x[i] + (i==theta ?delti[theta]:0);
|
| |
}
|
| |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
|
| |
|
| |
cptj=0;
|
| |
for(j=1; j<= nlstate; j++){
|
| |
for(i=1; i<=nlstate; i++){
|
| |
cptj=cptj+1;
|
| |
for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){
|
| |
gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
|
| |
}
|
| |
}
|
| |
}
|
| |
|
| |
|
| |
for(i=1; i<=npar; i++)
|
| |
xp[i] = x[i] - (i==theta ?delti[theta]:0);
|
| |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
|
| |
|
| |
cptj=0;
|
| |
for(j=1; j<= nlstate; j++){
|
| |
for(i=1;i<=nlstate;i++){
|
| |
cptj=cptj+1;
|
| |
for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){
|
| |
gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
|
| |
}
|
| |
}
|
| |
}
|
| |
|
| |
|
| |
|
| |
for(j=1; j<= nlstate*2; j++)
|
| |
for(h=0; h<=nhstepm-1; h++){
|
| |
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
|
| |
}
|
| |
|
| |
}
|
| |
|
| |
/* End theta */
|
| |
|
| |
trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);
|
| |
|
| |
for(h=0; h<=nhstepm-1; h++)
|
| |
for(j=1; j<=nlstate*2;j++)
|
| |
for(theta=1; theta <=npar; theta++)
|
| |
trgradg[h][j][theta]=gradg[h][theta][j];
|
| |
|
| |
|
| |
for(i=1;i<=nlstate*2;i++)
|
| |
for(j=1;j<=nlstate*2;j++)
|
| |
varhe[i][j][(int)age] =0.;
|
| |
|
| |
printf("%d|",(int)age);fflush(stdout);
|
| |
for(h=0;h<=nhstepm-1;h++){
|
| |
for(k=0;k<=nhstepm-1;k++){
|
| |
matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
|
| |
matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);
|
| |
for(i=1;i<=nlstate*2;i++)
|
| |
for(j=1;j<=nlstate*2;j++)
|
| |
varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
|
| |
}
|
| |
}
|
| |
|
| |
|
| |
/* Computing expectancies */
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| for(j=1; j<=nlstate;j++)
|
for(j=1; j<=nlstate;j++)
|
| for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
|
for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
|
| eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
|
eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
|
| /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
|
|
| |
/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
|
| |
|
| }
|
}
|
| |
|
| fprintf(ficreseij,"%3.0f",age );
|
fprintf(ficreseij,"%3.0f",age );
|
| |
cptj=0;
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| for(j=1; j<=nlstate;j++){
|
for(j=1; j<=nlstate;j++){
|
| fprintf(ficreseij," %9.4f", eij[i][j][(int)age]);
|
cptj++;
|
| |
fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );
|
| }
|
}
|
| fprintf(ficreseij,"\n");
|
fprintf(ficreseij,"\n");
|
| |
|
| |
free_matrix(gm,0,nhstepm,1,nlstate*2);
|
| |
free_matrix(gp,0,nhstepm,1,nlstate*2);
|
| |
free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);
|
| |
free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);
|
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
| }
|
}
|
| |
free_vector(xp,1,npar);
|
| |
free_matrix(dnewm,1,nlstate*2,1,npar);
|
| |
free_matrix(doldm,1,nlstate*2,1,nlstate*2);
|
| |
free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage);
|
| }
|
}
|
|
|
|
| /************ Variance ******************/
|
/************ Variance ******************/
|
| void varevsij(char fileres[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij)
|
void varevsij(char fileres[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm)
|
| {
|
{
|
| /* Variance of health expectancies */
|
/* Variance of health expectancies */
|
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
|
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
|
| double **newm;
|
double **newm;
|
| double **dnewm,**doldm;
|
double **dnewm,**doldm;
|
| int i, j, nhstepm, hstepm, h;
|
int i, j, nhstepm, hstepm, h, nstepm ;
|
| int k, cptcode;
|
int k, cptcode;
|
| double *xp;
|
double *xp;
|
| double **gp, **gm;
|
double **gp, **gm;
|
| double ***gradg, ***trgradg;
|
double ***gradg, ***trgradg;
|
| double ***p3mat;
|
double ***p3mat;
|
| double age,agelim;
|
double age,agelim, hf;
|
| int theta;
|
int theta;
|
|
|
|
| fprintf(ficresvij,"# Covariances of life expectancies\n");
|
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");
|
fprintf(ficresvij,"# Age");
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| for(j=1; j<=nlstate;j++)
|
for(j=1; j<=nlstate;j++)
|
|
Line 1609 void varevsij(char fileres[], double ***
|
Line 1735 void varevsij(char fileres[], double ***
|
| dnewm=matrix(1,nlstate,1,npar);
|
dnewm=matrix(1,nlstate,1,npar);
|
| doldm=matrix(1,nlstate,1,nlstate);
|
doldm=matrix(1,nlstate,1,nlstate);
|
|
|
|
| hstepm=1*YEARM; /* Every year of age */
|
if(estepm < stepm){
|
| hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */
|
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 only each 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;
|
agelim = AGESUP;
|
| for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
|
for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
|
| nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
|
nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
|
| if (stepm >= YEARM) hstepm=1;
|
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
|
| nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
|
|
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
| gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
|
gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
|
| gp=matrix(0,nhstepm,1,nlstate);
|
gp=matrix(0,nhstepm,1,nlstate);
|
|
Line 1670 void varevsij(char fileres[], double ***
|
Line 1810 void varevsij(char fileres[], double ***
|
| for(theta=1; theta <=npar; theta++)
|
for(theta=1; theta <=npar; theta++)
|
| trgradg[h][j][theta]=gradg[h][theta][j];
|
trgradg[h][j][theta]=gradg[h][theta][j];
|
|
|
|
| |
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
|
| for(i=1;i<=nlstate;i++)
|
for(i=1;i<=nlstate;i++)
|
| for(j=1;j<=nlstate;j++)
|
for(j=1;j<=nlstate;j++)
|
| vareij[i][j][(int)age] =0.;
|
vareij[i][j][(int)age] =0.;
|
| |
|
| for(h=0;h<=nhstepm;h++){
|
for(h=0;h<=nhstepm;h++){
|
| for(k=0;k<=nhstepm;k++){
|
for(k=0;k<=nhstepm;k++){
|
| matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
|
matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
|
| matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
|
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
|
| for(i=1;i<=nlstate;i++)
|
for(i=1;i<=nlstate;i++)
|
| for(j=1;j<=nlstate;j++)
|
for(j=1;j<=nlstate;j++)
|
| vareij[i][j][(int)age] += doldm[i][j];
|
vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
|
| }
|
}
|
| }
|
}
|
| h=1;
|
|
| if (stepm >= YEARM) h=stepm/YEARM;
|
|
| fprintf(ficresvij,"%.0f ",age );
|
fprintf(ficresvij,"%.0f ",age );
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| for(j=1; j<=nlstate;j++){
|
for(j=1; j<=nlstate;j++){
|
| fprintf(ficresvij," %.4f", h*vareij[i][j][(int)age]);
|
fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);
|
| }
|
}
|
| fprintf(ficresvij,"\n");
|
fprintf(ficresvij,"\n");
|
| free_matrix(gp,0,nhstepm,1,nlstate);
|
free_matrix(gp,0,nhstepm,1,nlstate);
|
|
Line 1718 void varprevlim(char fileres[], double *
|
Line 1859 void varprevlim(char fileres[], double *
|
| double age,agelim;
|
double age,agelim;
|
| int theta;
|
int theta;
|
|
|
|
| fprintf(ficresvpl,"# Standard deviation of prevalences limit\n");
|
fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n");
|
| fprintf(ficresvpl,"# Age");
|
fprintf(ficresvpl,"# Age");
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| fprintf(ficresvpl," %1d-%1d",i,i);
|
fprintf(ficresvpl," %1d-%1d",i,i);
|
|
Line 1787 void varprevlim(char fileres[], double *
|
Line 1928 void varprevlim(char fileres[], double *
|
| }
|
}
|
|
|
|
| /************ Variance of one-step probabilities ******************/
|
/************ Variance of one-step probabilities ******************/
|
| void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)
|
void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
|
| {
|
{
|
| int i, j;
|
int i, j, i1, k1, j1, z1;
|
| int k=0, cptcode;
|
int k=0, cptcode;
|
| double **dnewm,**doldm;
|
double **dnewm,**doldm;
|
| double *xp;
|
double *xp;
|
|
Line 1804 void varprob(char fileres[], double **ma
|
Line 1945 void varprob(char fileres[], double **ma
|
| if((ficresprob=fopen(fileresprob,"w"))==NULL) {
|
if((ficresprob=fopen(fileresprob,"w"))==NULL) {
|
| printf("Problem with resultfile: %s\n", fileresprob);
|
printf("Problem with resultfile: %s\n", fileresprob);
|
| }
|
}
|
| printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);
|
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
|
|
|
|
| |
fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");
|
| |
fprintf(ficresprob,"# Age");
|
| |
for(i=1; i<=nlstate;i++)
|
| |
for(j=1; j<=(nlstate+ndeath);j++)
|
| |
fprintf(ficresprob," p%1d-%1d (SE)",i,j);
|
| |
|
| |
|
| |
fprintf(ficresprob,"\n");
|
| |
|
|
|
|
| xp=vector(1,npar);
|
xp=vector(1,npar);
|
| dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
|
dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
|
| doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
|
doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
|
|
|
|
| cov[1]=1;
|
cov[1]=1;
|
| for (age=bage; age<=fage; age ++){
|
j=cptcoveff;
|
| cov[2]=age;
|
if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
| gradg=matrix(1,npar,1,9);
|
j1=0;
|
| trgradg=matrix(1,9,1,npar);
|
for(k1=1; k1<=1;k1++){
|
| gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
| gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
|
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#");
|
| |
}
|
|
|
|
| for(theta=1; theta <=npar; theta++){
|
for (age=bage; age<=fage; age ++){
|
| for(i=1; i<=npar; i++)
|
cov[2]=age;
|
| xp[i] = x[i] + (i==theta ?delti[theta]:0);
|
for (k=1; k<=cptcovn;k++) {
|
|
|
cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
|
| pmij(pmmij,cov,ncovmodel,xp,nlstate);
|
|
|
|
|
| k=0;
|
|
| for(i=1; i<= (nlstate+ndeath); i++){
|
|
| for(j=1; j<=(nlstate+ndeath);j++){
|
|
| k=k+1;
|
|
| gp[k]=pmmij[i][j];
|
|
| }
|
}
|
| }
|
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
|
|
|
for (k=1; k<=cptcovprod;k++)
|
| for(i=1; i<=npar; i++)
|
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
|
| xp[i] = x[i] - (i==theta ?delti[theta]:0);
|
|
| |
gradg=matrix(1,npar,1,9);
|
| |
trgradg=matrix(1,9,1,npar);
|
| |
gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
|
| |
gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
|
| |
|
| |
for(theta=1; theta <=npar; theta++){
|
| |
for(i=1; i<=npar; i++)
|
| |
xp[i] = x[i] + (i==theta ?delti[theta]:0);
|
| |
|
| |
pmij(pmmij,cov,ncovmodel,xp,nlstate);
|
| |
|
| |
k=0;
|
| |
for(i=1; i<= (nlstate+ndeath); 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]:0);
|
|
|
|
|
|
pmij(pmmij,cov,ncovmodel,xp,nlstate);
|
| pmij(pmmij,cov,ncovmodel,xp,nlstate);
|
k=0;
|
| k=0;
|
for(i=1; i<=(nlstate+ndeath); i++){
|
| for(i=1; i<=(nlstate+ndeath); i++){
|
for(j=1; j<=(nlstate+ndeath);j++){
|
| for(j=1; j<=(nlstate+ndeath);j++){
|
k=k+1;
|
| k=k+1;
|
gm[k]=pmmij[i][j];
|
| gm[k]=pmmij[i][j];
|
}
|
| }
|
}
|
| }
|
|
|
|
|
| for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
|
for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
|
| gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];
|
gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];
|
| }
|
}
|
|
|
|
| for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
|
|
| for(theta=1; theta <=npar; theta++)
|
|
| trgradg[j][theta]=gradg[theta][j];
|
|
|
|
|
| matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
|
|
| matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
|
|
|
|
|
| pmij(pmmij,cov,ncovmodel,x,nlstate);
|
|
|
|
|
| k=0;
|
for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
|
| for(i=1; i<=(nlstate+ndeath); i++){
|
for(theta=1; theta <=npar; theta++)
|
| for(j=1; j<=(nlstate+ndeath);j++){
|
trgradg[j][theta]=gradg[theta][j];
|
| k=k+1;
|
|
| gm[k]=pmmij[i][j];
|
matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
|
| |
matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
|
| |
|
| |
pmij(pmmij,cov,ncovmodel,x,nlstate);
|
| |
|
| |
k=0;
|
| |
for(i=1; i<=(nlstate+ndeath); i++){
|
| |
for(j=1; j<=(nlstate+ndeath);j++){
|
| |
k=k+1;
|
| |
gm[k]=pmmij[i][j];
|
| |
}
|
| }
|
}
|
| }
|
|
|
|
|
| /*printf("\n%d ",(int)age);
|
/*printf("\n%d ",(int)age);
|
| for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
|
for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
|
|
|
|
|
|
|
| printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
|
printf("%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(ficresprob,"\n%d ",(int)age);
|
|
|
|
| for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
|
|
| if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
|
|
| if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
|
|
| }
|
|
|
|
|
| |
for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
|
| |
fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));
|
| |
|
| |
}
|
| |
}
|
| free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
|
free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
|
| free_vector(gm,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(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
|
| free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
|
free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
|
| }
|
}
|
| free_vector(xp,1,npar);
|
free_vector(xp,1,npar);
|
| fclose(ficresprob);
|
fclose(ficresprob);
|
|
|
|
| }
|
}
|
|
|
|
| /******************* Printing html file ***********/
|
/******************* Printing html file ***********/
|
| void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[],char rfileres[] ){
|
void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
|
| |
int lastpass, int stepm, int weightopt, char model[],\
|
| |
int imx,int jmin, int jmax, double jmeanint,char optionfile[], \
|
| |
char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\
|
| |
char version[], int popforecast, int estepm ,/* \ */
|
| |
double jprev1, double mprev1,double anprev1, \
|
| |
double jprev2, double mprev2,double anprev2){
|
| int jj1, k1, i1, cpt;
|
int jj1, k1, i1, cpt;
|
| FILE *fichtm;
|
FILE *fichtm;
|
| /*char optionfilehtm[FILENAMELENGTH];*/
|
/*char optionfilehtm[FILENAMELENGTH];*/
|
|
Line 1903 void printinghtml(char fileres[], char t
|
Line 2077 void printinghtml(char fileres[], char t
|
| printf("Problem with %s \n",optionfilehtm), exit(0);
|
printf("Problem with %s \n",optionfilehtm), exit(0);
|
| }
|
}
|
|
|
|
| fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.8 </font> <hr size=\"2\" color=\"#EC5E5E\">
|
fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n
|
| Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
|
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
|
|
|
\n
|
| Total number of observations=%d <br>
|
Total number of observations=%d <br>\n
|
| Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
|
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
|
| <hr size=\"2\" color=\"#EC5E5E\">
|
<hr size=\"2\" color=\"#EC5E5E\">
|
| <li>Outputs files<br><br>\n
|
<ul><li>Parameter files<br>\n
|
| - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
|
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
|
| - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>
|
- Gnuplot file name: <a href=\"%s\">%s</a><br></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot);
|
| - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>
|
|
| - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>
|
fprintf(fichtm,"<ul><li>Result files (first order: no variance)<br>\n
|
| - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>
|
- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n
|
| - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
|
- Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n
|
| - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
|
- Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n
|
| - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
|
- Life expectancies by age and initial health status (estepm=%2d months):
|
| - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
|
<a href=\"e%s\">e%s</a> <br>\n</li>", \
|
| - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>
|
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
|
| - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>
|
|
| <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
|
fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n
|
|
|
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n
|
| |
- Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
|
| |
- Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
|
| |
- Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
|
| |
- Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
|
| |
|
| |
if(popforecast==1) fprintf(fichtm,"\n
|
| |
- Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
|
| |
- Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n
|
| |
<br>",fileres,fileres,fileres,fileres);
|
| |
else
|
| |
fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);
|
| fprintf(fichtm," <li>Graphs</li><p>");
|
fprintf(fichtm," <li>Graphs</li><p>");
|
|
|
|
| m=cptcoveff;
|
m=cptcoveff;
|
|
Line 1938 fprintf(fichtm," <li>Graphs</li><p>");
|
Line 2123 fprintf(fichtm," <li>Graphs</li><p>");
|
| fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
|
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
|
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
|
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
|
| }
|
}
|
| fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
|
/* Pij */
|
| <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
|
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>
|
| |
<img src=\"pe%s%d1.png\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);
|
| |
/* Quasi-incidences */
|
| |
fprintf(fichtm,"<br>- 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: pe%s%d2.png<br>
|
| |
<img src=\"pe%s%d2.png\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);
|
| |
/* Stable prevalence in each health state */
|
| for(cpt=1; cpt<nlstate;cpt++){
|
for(cpt=1; cpt<nlstate;cpt++){
|
| fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>
|
fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
|
| <img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
| }
|
}
|
| for(cpt=1; cpt<=nlstate;cpt++) {
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
| fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
|
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
|
| interval) in state (%d): v%s%d%d.gif <br>
|
interval) in state (%d): v%s%d%d.png <br>
|
| <img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
| }
|
}
|
| for(cpt=1; cpt<=nlstate;cpt++) {
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
| fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>
|
fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
|
| <img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
| }
|
}
|
| fprintf(fichtm,"\n<br>- Total life expectancy by age and
|
fprintf(fichtm,"\n<br>- Total life expectancy by age and
|
| health expectancies in states (1) and (2): e%s%d.gif<br>
|
health expectancies in states (1) and (2): e%s%d.png<br>
|
| <img src=\"e%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
|
<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
|
| fprintf(fichtm,"\n</body>");
|
fprintf(fichtm,"\n</body>");
|
| }
|
}
|
| }
|
}
|
|
Line 1963 fclose(fichtm);
|
Line 2153 fclose(fichtm);
|
| }
|
}
|
|
|
|
| /******************* Gnuplot file **************/
|
/******************* Gnuplot file **************/
|
| void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemaxpar, double fage , char pathc[], double p[]){
|
void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
|
|
|
|
| int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
|
int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
|
|
|
int ng;
|
| strcpy(optionfilegnuplot,optionfilefiname);
|
strcpy(optionfilegnuplot,optionfilefiname);
|
| strcat(optionfilegnuplot,".plt");
|
strcat(optionfilegnuplot,".gp.txt");
|
| if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
|
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
|
| printf("Problem with file %s",optionfilegnuplot);
|
printf("Problem with file %s",optionfilegnuplot);
|
| }
|
}
|
|
Line 1983 m=pow(2,cptcoveff);
|
Line 2173 m=pow(2,cptcoveff);
|
| for (k1=1; k1<= m ; k1 ++) {
|
for (k1=1; k1<= m ; k1 ++) {
|
|
|
|
| #ifdef windows
|
#ifdef windows
|
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",agemin,fage,fileres,k1-1,k1-1);
|
fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
|
| |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
|
| #endif
|
#endif
|
| #ifdef unix
|
#ifdef unix
|
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres);
|
fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
|
| |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
|
| #endif
|
#endif
|
|
|
|
| for (i=1; i<= nlstate ; i ++) {
|
for (i=1; i<= nlstate ; i ++) {
|
|
Line 2005 for (i=1; i<= nlstate ; i ++) {
|
Line 2197 for (i=1; i<= nlstate ; i ++) {
|
| }
|
}
|
| fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
|
fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
|
| #ifdef unix
|
#ifdef unix
|
| fprintf(ficgp,"\nset ter gif small size 400,300");
|
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n");
|
| #endif
|
#endif
|
| fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
| }
|
}
|
| }
|
}
|
| /*2 eme*/
|
/*2 eme*/
|
|
|
|
| for (k1=1; k1<= m ; k1 ++) {
|
for (k1=1; k1<= m ; k1 ++) {
|
| fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);
|
fprintf(ficgp,"\nset out \"e%s%d.png\" \n",strtok(optionfile, "."),k1);
|
| |
fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
|
|
|
|
| for (i=1; i<= nlstate+1 ; i ++) {
|
for (i=1; i<= nlstate+1 ; i ++) {
|
| k=2*i;
|
k=2*i;
|
|
Line 2038 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
|
Line 2230 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
|
| if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
|
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
|
| else fprintf(ficgp,"\" t\"\" w l 0,");
|
else fprintf(ficgp,"\" t\"\" w l 0,");
|
| }
|
}
|
| fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);
|
|
| }
|
}
|
|
|
|
| /*3eme*/
|
/*3eme*/
|
|
|
|
| for (k1=1; k1<= m ; k1 ++) {
|
for (k1=1; k1<= m ; k1 ++) {
|
| for (cpt=1; cpt<= nlstate ; cpt ++) {
|
for (cpt=1; cpt<= nlstate ; cpt ++) {
|
| k=2+nlstate*(cpt-1);
|
k=2+nlstate*(2*cpt-2);
|
| fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k1-1,k1-1,k,cpt);
|
fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
|
| |
fprintf(ficgp,"set ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
|
| |
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
|
| |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
|
| |
fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
|
| |
fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
|
| |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
|
| |
fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
|
| |
|
| |
*/
|
| for (i=1; i< nlstate ; i ++) {
|
for (i=1; i< nlstate ; i ++) {
|
| fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);
|
fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1);
|
| |
|
| }
|
}
|
| fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
| }
|
|
| }
|
}
|
| |
}
|
|
|
|
| /* CV preval stat */
|
/* CV preval stat */
|
| for (k1=1; k1<= m ; k1 ++) {
|
for (k1=1; k1<= m ; k1 ++) {
|
| for (cpt=1; cpt<nlstate ; cpt ++) {
|
for (cpt=1; cpt<nlstate ; cpt ++) {
|
| k=3;
|
k=3;
|
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemaxpar,fileres,k1,k+cpt+1,k+1);
|
fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
|
| |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
|
|
|
|
| for (i=1; i< nlstate ; i ++)
|
for (i=1; i< nlstate ; i ++)
|
| fprintf(ficgp,"+$%d",k+i+1);
|
fprintf(ficgp,"+$%d",k+i+1);
|
|
Line 2071 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
|
Line 2272 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
|
| fprintf(ficgp,"+$%d",l+i+1);
|
fprintf(ficgp,"+$%d",l+i+1);
|
| }
|
}
|
| fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);
|
fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);
|
| fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
| }
|
}
|
| }
|
}
|
|
|
|
|
Line 2087 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
|
Line 2287 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
|
|
|
| for(jk=1; jk <=m; jk++) {
|
for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
|
| fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot [%.f:%.f] ",agemin,agemaxpar);
|
for(jk=1; jk <=m; jk++) {
|
| i=1;
|
fprintf(ficgp,"\nset out \"pe%s%d%d.png\" \n",strtok(optionfile, "."),jk,ng);
|
| for(k2=1; k2<=nlstate; k2++) {
|
if (ng==2)
|
| k3=i;
|
fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
|
| for(k=1; k<=(nlstate+ndeath); k++) {
|
else
|
| if (k != k2){
|
fprintf(ficgp,"\nset title \"Probability\"\n");
|
| fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
|
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar);
|
| ij=1;
|
i=1;
|
| for(j=3; j <=ncovmodel; j++) {
|
for(k2=1; k2<=nlstate; k2++) {
|
| if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
k3=i;
|
| fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
for(k=1; k<=(nlstate+ndeath); k++) {
|
| ij++;
|
if (k != k2){
|
| }
|
if(ng==2)
|
| else
|
fprintf(ficgp," %f*exp(p%d+p%d*x",stepm/YEARM,i,i+1);
|
| fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
else
|
| }
|
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
|
| fprintf(ficgp,")/(1");
|
ij=1;
|
|
|
for(j=3; j <=ncovmodel; j++) {
|
| for(k1=1; k1 <=nlstate; k1++){
|
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
| fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
|
fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
| ij=1;
|
ij++;
|
| for(j=3; j <=ncovmodel; j++){
|
}
|
| if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
else
|
| fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
| ij++;
|
}
|
| }
|
fprintf(ficgp,")/(1");
|
| else
|
|
| fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
for(k1=1; k1 <=nlstate; k1++){
|
| }
|
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
|
| fprintf(ficgp,")");
|
ij=1;
|
| }
|
for(j=3; j <=ncovmodel; j++){
|
| fprintf(ficgp,") t \"p%d%d\" ", k2,k);
|
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
| if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
|
fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
| i=i+ncovmodel;
|
ij++;
|
| |
}
|
| |
else
|
| |
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
| |
}
|
| |
fprintf(ficgp,")");
|
| |
}
|
| |
fprintf(ficgp,") t \"p%d%d\" ", k2,k);
|
| |
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
|
| |
i=i+ncovmodel;
|
| |
}
|
| |
}
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
| fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
|
fclose(ficgp);
|
| }
|
|
|
|
|
| fclose(ficgp);
|
|
| } /* end gnuplot */
|
} /* end gnuplot */
|
|
|
|
|
|
|
| /*************** Moving average **************/
|
/*************** Moving average **************/
|
| void movingaverage(double agedeb, double fage,double agemin, double ***mobaverage){
|
void movingaverage(double agedeb, double fage,double ageminpar, double ***mobaverage){
|
|
|
|
| int i, cpt, cptcod;
|
int i, cpt, cptcod;
|
| for (agedeb=agemin; agedeb<=fage; agedeb++)
|
for (agedeb=ageminpar; agedeb<=fage; agedeb++)
|
| for (i=1; i<=nlstate;i++)
|
for (i=1; i<=nlstate;i++)
|
| for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
|
for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
|
| mobaverage[(int)agedeb][i][cptcod]=0.;
|
mobaverage[(int)agedeb][i][cptcod]=0.;
|
|
|
|
| for (agedeb=agemin+4; agedeb<=fage; agedeb++){
|
for (agedeb=ageminpar+4; agedeb<=fage; agedeb++){
|
| for (i=1; i<=nlstate;i++){
|
for (i=1; i<=nlstate;i++){
|
| for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
|
for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
|
| for (cpt=0;cpt<=4;cpt++){
|
for (cpt=0;cpt<=4;cpt++){
|
|
Line 2158 void movingaverage(double agedeb, double
|
Line 2366 void movingaverage(double agedeb, double
|
|
|
|
|
|
|
| /************** Forecasting ******************/
|
/************** Forecasting ******************/
|
| prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){
|
prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){
|
|
|
|
| int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
|
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
|
| int *popage;
|
int *popage;
|
|
Line 2170 prevforecast(char fileres[], double anpr
|
Line 2378 prevforecast(char fileres[], double anpr
|
| agelim=AGESUP;
|
agelim=AGESUP;
|
| calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
|
calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
|
|
|
|
| prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
|
|
|
|
|
|
| strcpy(fileresf,"f");
|
strcpy(fileresf,"f");
|
|
Line 2184 calagedate=(anproj1+mproj1/12.+jproj1/36
|
Line 2392 calagedate=(anproj1+mproj1/12.+jproj1/36
|
|
|
|
| if (mobilav==1) {
|
if (mobilav==1) {
|
| mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
| movingaverage(agedeb, fage, agemin, mobaverage);
|
movingaverage(agedeb, fage, ageminpar, mobaverage);
|
| }
|
}
|
|
|
|
| stepsize=(int) (stepm+YEARM-1)/YEARM;
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
|
Line 2221 calagedate=(anproj1+mproj1/12.+jproj1/36
|
Line 2429 calagedate=(anproj1+mproj1/12.+jproj1/36
|
| fprintf(ficresf,"\n");
|
fprintf(ficresf,"\n");
|
| fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);
|
fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);
|
|
|
|
| for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
|
for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){
|
| nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
| nhstepm = nhstepm/hstepm;
|
nhstepm = nhstepm/hstepm;
|
|
|
|
|
Line 2231 calagedate=(anproj1+mproj1/12.+jproj1/36
|
Line 2439 calagedate=(anproj1+mproj1/12.+jproj1/36
|
|
|
|
| for (h=0; h<=nhstepm; h++){
|
for (h=0; h<=nhstepm; h++){
|
| if (h==(int) (calagedate+YEARM*cpt)) {
|
if (h==(int) (calagedate+YEARM*cpt)) {
|
| fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
|
fprintf(ficresf,"\n %.f %.f ",anproj1+cpt,agedeb+h*hstepm/YEARM*stepm);
|
| }
|
}
|
| for(j=1; j<=nlstate+ndeath;j++) {
|
for(j=1; j<=nlstate+ndeath;j++) {
|
| kk1=0.;kk2=0;
|
kk1=0.;kk2=0;
|
|
Line 2260 calagedate=(anproj1+mproj1/12.+jproj1/36
|
Line 2468 calagedate=(anproj1+mproj1/12.+jproj1/36
|
| fclose(ficresf);
|
fclose(ficresf);
|
| }
|
}
|
| /************** Forecasting ******************/
|
/************** Forecasting ******************/
|
| populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
|
populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
|
|
|
|
| int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
|
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
|
| int *popage;
|
int *popage;
|
|
Line 2274 populforecast(char fileres[], double anp
|
Line 2482 populforecast(char fileres[], double anp
|
| agelim=AGESUP;
|
agelim=AGESUP;
|
| calagedate=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
|
calagedate=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
|
|
|
|
| prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
|
|
|
|
|
|
| strcpy(filerespop,"pop");
|
strcpy(filerespop,"pop");
|
|
Line 2288 populforecast(char fileres[], double anp
|
Line 2496 populforecast(char fileres[], double anp
|
|
|
|
| if (mobilav==1) {
|
if (mobilav==1) {
|
| mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
| movingaverage(agedeb, fage, agemin, mobaverage);
|
movingaverage(agedeb, fage, ageminpar, mobaverage);
|
| }
|
}
|
|
|
|
| stepsize=(int) (stepm+YEARM-1)/YEARM;
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
|
Line 2329 populforecast(char fileres[], double anp
|
Line 2537 populforecast(char fileres[], double anp
|
| for (cpt=0; cpt<=0;cpt++) {
|
for (cpt=0; cpt<=0;cpt++) {
|
| fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
|
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
|
|
|
|
| for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
|
for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){
|
| nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
| nhstepm = nhstepm/hstepm;
|
nhstepm = nhstepm/hstepm;
|
|
|
|
|
Line 2375 populforecast(char fileres[], double anp
|
Line 2583 populforecast(char fileres[], double anp
|
|
|
|
| for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
|
for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
|
| fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
|
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
|
| for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
|
for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){
|
| nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
| nhstepm = nhstepm/hstepm;
|
nhstepm = nhstepm/hstepm;
|
|
|
|
|
Line 2421 int main(int argc, char *argv[])
|
Line 2629 int main(int argc, char *argv[])
|
|
|
|
| int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
|
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
|
| double agedeb, agefin,hf;
|
double agedeb, agefin,hf;
|
| double agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
|
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
|
|
|
|
| double fret;
|
double fret;
|
| double **xi,tmp,delta;
|
double **xi,tmp,delta;
|
|
Line 2448 int main(int argc, char *argv[])
|
Line 2656 int main(int argc, char *argv[])
|
| int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
|
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
|
| int mobilav=0,popforecast=0;
|
int mobilav=0,popforecast=0;
|
| int hstepm, nhstepm;
|
int hstepm, nhstepm;
|
| double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;
|
double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;
|
|
|
|
| double bage, fage, age, agelim, agebase;
|
double bage, fage, age, agelim, agebase;
|
| double ftolpl=FTOL;
|
double ftolpl=FTOL;
|
|
Line 2466 int main(int argc, char *argv[])
|
Line 2674 int main(int argc, char *argv[])
|
| double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
|
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
|
|
|
|
|
|
|
| char version[80]="Imach version 0.8, March 2002, INED-EUROREVES ";
|
char version[80]="Imach version 0.8e, May 2002, INED-EUROREVES ";
|
| char *alph[]={"a","a","b","c","d","e"}, str[4];
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
|
|
|
|
|
|
|
Line 2479 int main(int argc, char *argv[])
|
Line 2687 int main(int argc, char *argv[])
|
| struct timeval start_time, end_time;
|
struct timeval start_time, end_time;
|
|
|
|
| gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
|
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
|
|
|
getcwd(pathcd, size);
|
|
|
|
| printf("\n%s",version);
|
printf("\n%s",version);
|
| if(argc <=1){
|
if(argc <=1){
|
|
Line 2705 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 2913 while((c=getc(ficpar))=='#' && c!= EOF){
|
| if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
|
if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
|
| if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
|
if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
|
| if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
|
if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
|
| }
|
}*/
|
|
|
/* for (i=1; i<=imx; i++){
|
| for (i=1; i<=imx; i++)
|
if (s[4][i]==9) s[4][i]=-1;
|
| if (covar[1][i]==0) printf("%d %.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]));*/
|
printf("%d %.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]));}*/
|
|
|
|
| |
|
| /* Calculation of the number of parameter from char model*/
|
/* Calculation of the number of parameter from char model*/
|
| Tvar=ivector(1,15);
|
Tvar=ivector(1,15);
|
| Tprod=ivector(1,15);
|
Tprod=ivector(1,15);
|
|
Line 2724 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 2933 while((c=getc(ficpar))=='#' && c!= EOF){
|
| cptcovn=j+1;
|
cptcovn=j+1;
|
| cptcovprod=j1;
|
cptcovprod=j1;
|
|
|
|
|
|
|
| strcpy(modelsav,model);
|
strcpy(modelsav,model);
|
| if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
|
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
|
| printf("Error. Non available option model=%s ",model);
|
printf("Error. Non available option model=%s ",model);
|
|
Line 2780 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 2988 while((c=getc(ficpar))=='#' && c!= EOF){
|
| }
|
}
|
| }
|
}
|
|
|
|
| /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
|
/* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
|
| printf("cptcovprod=%d ", cptcovprod);
|
printf("cptcovprod=%d ", cptcovprod);
|
| scanf("%d ",i);*/
|
scanf("%d ",i);*/
|
| fclose(fic);
|
fclose(fic);
|
|
Line 2792 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3000 while((c=getc(ficpar))=='#' && c!= EOF){
|
| /*-calculation of age at interview from date of interview and age at death -*/
|
/*-calculation of age at interview from date of interview and age at death -*/
|
| agev=matrix(1,maxwav,1,imx);
|
agev=matrix(1,maxwav,1,imx);
|
|
|
|
| for (i=1; i<=imx; i++)
|
for (i=1; i<=imx; i++) {
|
| for(m=2; (m<= maxwav); m++)
|
for(m=2; (m<= maxwav); m++) {
|
| if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
|
if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
|
| anint[m][i]=9999;
|
anint[m][i]=9999;
|
| s[m][i]=-1;
|
s[m][i]=-1;
|
| }
|
}
|
|
|
if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;
|
| |
}
|
| |
}
|
| |
|
| for (i=1; i<=imx; i++) {
|
for (i=1; i<=imx; i++) {
|
| agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
|
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
|
| for(m=1; (m<= maxwav); m++){
|
for(m=1; (m<= maxwav); m++){
|
| if(s[m][i] >0){
|
if(s[m][i] >0){
|
| if (s[m][i] == nlstate+1) {
|
if (s[m][i] >= nlstate+1) {
|
| if(agedc[i]>0)
|
if(agedc[i]>0)
|
| if(moisdc[i]!=99 && andc[i]!=9999)
|
if(moisdc[i]!=99 && andc[i]!=9999)
|
| agev[m][i]=agedc[i];
|
agev[m][i]=agedc[i];
|
| else {
|
/*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
|
| |
else {
|
| if (andc[i]!=9999){
|
if (andc[i]!=9999){
|
| printf("Warning negative age at death: %d line:%d\n",num[i],i);
|
printf("Warning negative age at death: %d line:%d\n",num[i],i);
|
| agev[m][i]=-1;
|
agev[m][i]=-1;
|
|
Line 2882 printf("Total number of individuals= %d,
|
Line 3094 printf("Total number of individuals= %d,
|
| for(j=1; j <= ncodemax[k]; j++){
|
for(j=1; j <= ncodemax[k]; j++){
|
| for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
|
for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
|
| h++;
|
h++;
|
| if (h>m) h=1;codtab[h][k]=j;
|
if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
|
| |
/* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
|
|
/* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]);
|
|
|
codtab[1][2]=1;codtab[2][2]=2; */
|
| /*for(i=1; i <=m ;i++){
|
/* for(i=1; i <=m ;i++){
|
| for(k=1; k <=cptcovn; k++){
|
for(k=1; k <=cptcovn; k++){
|
| printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);
|
printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
|
| }
|
}
|
| printf("\n");
|
printf("\n");
|
| }
|
}
|
| scanf("%d",i);*/
|
scanf("%d",i);*/
|
|
|
|
| /* Calculates basic frequencies. Computes observed prevalence at single age
|
/* Calculates basic frequencies. Computes observed prevalence at single age
|
| and prints on file fileres'p'. */
|
and prints on file fileres'p'. */
|
|
Line 2988 printf("Total number of individuals= %d,
|
Line 3201 printf("Total number of individuals= %d,
|
| fputs(line,ficparo);
|
fputs(line,ficparo);
|
| }
|
}
|
| ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
|
|
estepm=0;
|
| fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemaxpar, &bage, &fage);
|
fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
|
|
|
if (estepm==0 || estepm < stepm) estepm=stepm;
|
| if (fage <= 2) {
|
if (fage <= 2) {
|
| bage = agemin;
|
bage = ageminpar;
|
| fage = agemaxpar;
|
fage = agemaxpar;
|
| }
|
}
|
|
|
|
| fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
|
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
|
| fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);
|
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
|
| fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);
|
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
|
|
|
|
| while((c=getc(ficpar))=='#' && c!= EOF){
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
| ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
|
Line 3056 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3269 while((c=getc(ficpar))=='#' && c!= EOF){
|
| freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
|
|
|
| /*------------ gnuplot -------------*/
|
/*------------ gnuplot -------------*/
|
| printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemaxpar,fage, pathc,p);
|
printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, ageminpar,agemaxpar,fage, pathc,p);
|
|
|
|
| /*------------ free_vector -------------*/
|
/*------------ free_vector -------------*/
|
| chdir(path);
|
chdir(path);
|
|
Line 3072 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3285 while((c=getc(ficpar))=='#' && c!= EOF){
|
|
|
|
| /*--------- index.htm --------*/
|
/*--------- index.htm --------*/
|
|
|
|
| printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres);
|
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
|
|
|
|
|
|
|
| /*--------------- Prevalence limit --------------*/
|
/*--------------- Prevalence limit --------------*/
|
|
Line 3095 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3308 while((c=getc(ficpar))=='#' && c!= EOF){
|
| savms= 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 */
|
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
|
| k=0;
|
k=0;
|
| agebase=agemin;
|
agebase=ageminpar;
|
| agelim=agemaxpar;
|
agelim=agemaxpar;
|
| ftolpl=1.e-10;
|
ftolpl=1.e-10;
|
| i1=cptcoveff;
|
i1=cptcoveff;
|
|
Line 3156 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3369 while((c=getc(ficpar))=='#' && c!= EOF){
|
| for(j=1; j<=nlstate+ndeath;j++)
|
for(j=1; j<=nlstate+ndeath;j++)
|
| fprintf(ficrespij," %1d-%1d",i,j);
|
fprintf(ficrespij," %1d-%1d",i,j);
|
| fprintf(ficrespij,"\n");
|
fprintf(ficrespij,"\n");
|
| for (h=0; h<=nhstepm; h++){
|
for (h=0; h<=nhstepm; h++){
|
| fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
|
fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
|
| for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
| for(j=1; j<=nlstate+ndeath;j++)
|
for(j=1; j<=nlstate+ndeath;j++)
|
| fprintf(ficrespij," %.5f", p3mat[i][j][h]);
|
fprintf(ficrespij," %.5f", p3mat[i][j][h]);
|
| fprintf(ficrespij,"\n");
|
fprintf(ficrespij,"\n");
|
| }
|
}
|
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
| fprintf(ficrespij,"\n");
|
fprintf(ficrespij,"\n");
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
|
|
|
| /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
|
varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
|
|
|
|
| fclose(ficrespij);
|
fclose(ficrespij);
|
|
|
|
|
Line 3178 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3391 while((c=getc(ficpar))=='#' && c!= EOF){
|
| if((stepm == 1) && (strcmp(model,".")==0)){
|
if((stepm == 1) && (strcmp(model,".")==0)){
|
| prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
|
prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
|
| if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
|
if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
|
| free_matrix(mint,1,maxwav,1,n);
|
}
|
| free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
|
|
| free_vector(weight,1,n);}
|
|
| else{
|
else{
|
| erreur=108;
|
erreur=108;
|
| printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);
|
printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);
|
|
Line 3210 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3421 while((c=getc(ficpar))=='#' && c!= EOF){
|
| printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
|
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
|
| }
|
}
|
| printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
|
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
|
| |
calagedate=-1;
|
| |
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
|
|
|
| k=0;
|
k=0;
|
| for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
|
Line 3222 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3435 while((c=getc(ficpar))=='#' && c!= EOF){
|
|
|
|
| fprintf(ficreseij,"\n#****** ");
|
fprintf(ficreseij,"\n#****** ");
|
| for(j=1;j<=cptcoveff;j++)
|
for(j=1;j<=cptcoveff;j++)
|
| fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
|
fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
| fprintf(ficreseij,"******\n");
|
fprintf(ficreseij,"******\n");
|
|
|
|
| fprintf(ficresvij,"\n#****** ");
|
fprintf(ficresvij,"\n#****** ");
|
| for(j=1;j<=cptcoveff;j++)
|
for(j=1;j<=cptcoveff;j++)
|
| fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
|
fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
| fprintf(ficresvij,"******\n");
|
fprintf(ficresvij,"******\n");
|
|
|
|
| eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
|
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
|
| oldm=oldms;savm=savms;
|
oldm=oldms;savm=savms;
|
| evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k);
|
evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);
|
| |
|
| vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
|
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
|
| oldm=oldms;savm=savms;
|
oldm=oldms;savm=savms;
|
| varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
|
varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);
|
|
|
|
|
|
|
|
|
|
|
Line 3255 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3469 while((c=getc(ficpar))=='#' && c!= EOF){
|
| for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
|
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
|
| for(i=1, epj[j]=0.;i <=nlstate;i++) {
|
for(i=1, epj[j]=0.;i <=nlstate;i++) {
|
| epj[j] += prlim[i][i]*eij[i][j][(int)age];
|
epj[j] += prlim[i][i]*eij[i][j][(int)age];
|
| |
/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
|
| }
|
}
|
| epj[nlstate+1] +=epj[j];
|
epj[nlstate+1] +=epj[j];
|
| }
|
}
|
| |
|
| for(i=1, vepp=0.;i <=nlstate;i++)
|
for(i=1, vepp=0.;i <=nlstate;i++)
|
| for(j=1;j <=nlstate;j++)
|
for(j=1;j <=nlstate;j++)
|
| vepp += vareij[i][j][(int)age];
|
vepp += vareij[i][j][(int)age];
|
| fprintf(ficrest," %7.2f (%7.2f)", epj[nlstate+1],sqrt(vepp));
|
fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
|
| for(j=1;j <=nlstate;j++){
|
for(j=1;j <=nlstate;j++){
|
| fprintf(ficrest," %7.2f (%7.2f)", epj[j],sqrt(vareij[j][j][(int)age]));
|
fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
|
| }
|
}
|
| fprintf(ficrest,"\n");
|
fprintf(ficrest,"\n");
|
| }
|
}
|
| }
|
}
|
| }
|
}
|
|
|
free_matrix(mint,1,maxwav,1,n);
|
| |
free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
|
| |
free_vector(weight,1,n);
|
| fclose(ficreseij);
|
fclose(ficreseij);
|
| fclose(ficresvij);
|
fclose(ficresvij);
|
| fclose(ficrest);
|
fclose(ficrest);
|
|
Line 3345 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3563 while((c=getc(ficpar))=='#' && c!= EOF){
|
|
|
|
| #ifdef windows
|
#ifdef windows
|
| while (z[0] != 'q') {
|
while (z[0] != 'q') {
|
| chdir(path);
|
/* chdir(path); */
|
| printf("\nType e to edit output files, c to start again, and q for exiting: ");
|
printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
|
| scanf("%s",z);
|
scanf("%s",z);
|
| if (z[0] == 'c') system("./imach");
|
if (z[0] == 'c') system("./imach");
|
| else if (z[0] == 'e') {
|
else if (z[0] == 'e') system(optionfilehtm);
|
| chdir(path);
|
else if (z[0] == 'g') system(plotcmd);
|
| system(optionfilehtm);
|
|
| }
|
|
| else if (z[0] == 'q') exit(0);
|
else if (z[0] == 'q') exit(0);
|
| }
|
}
|
| #endif
|
#endif
|
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src