--- imach/src/imach.c 2002/02/20 17:02:08 1.13
+++ imach/src/imach.c 2002/02/20 17:08:52 1.15
@@ -73,6 +73,7 @@ int npar=NPARMAX;
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
int ncovmodel, ncov; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
+int popbased=0;
int *wav; /* Number of waves for this individuual 0 is possible */
int maxwav; /* Maxim number of waves */
@@ -1149,7 +1150,7 @@ void lubksb(double **a, int n, int *indx
}
/************ Frequencies ********************/
-void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax)
+void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1)
{ /* Some frequencies */
int i, m, jk, k1, i1, j1, bool, z1,z2,j;
@@ -1191,7 +1192,7 @@ void freqsummary(char fileres[], int ag
bool=0;
}
if (bool==1) {
- for(m=firstpass; m<=lastpass-1; m++){
+ for(m=fprev1; m<=lprev1; m++){
if(agev[m][i]==0) agev[m][i]=agemax+1;
if(agev[m][i]==1) agev[m][i]=agemax+2;
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
@@ -1215,7 +1216,7 @@ void freqsummary(char fileres[], int ag
printf("Age %d", i);
for(jk=1; jk <=nlstate ; jk++){
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(m=-1, pos=0; m <=0 ; m++)
@@ -1225,10 +1226,12 @@ void freqsummary(char fileres[], int ag
else
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
}
- for(jk=1; jk <=nlstate ; jk++){
- for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)
+
+ 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; jk <=nlstate ; jk++)
pos += pp[jk];
for(jk=1; jk <=nlstate ; jk++){
@@ -1262,6 +1265,83 @@ void freqsummary(char fileres[], int ag
} /* End of Freq */
+/************ Prevalence ********************/
+void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1)
+{ /* Some frequencies */
+
+ int i, m, jk, k1, i1, j1, bool, z1,z2,j;
+ double ***freq; /* Frequencies */
+ double *pp;
+ double pos;
+
+ pp=vector(1,nlstate);
+ probs= ma3x(1,130 ,1,8, 1,8);
+
+ freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
+ j1=0;
+
+ 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+ndeath; i++)
+ for (jk=-1; jk<=nlstate+ndeath; jk++)
+ for(m=agemin; m <= agemax+3; m++)
+ freq[i][jk][m]=0;
+
+ for (i=1; i<=imx; i++) {
+ bool=1;
+ if (cptcovn>0) {
+ for (z1=1; z1<=cptcoveff; z1++)
+ if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
+ bool=0;
+ }
+ if (bool==1) {
+ for(m=fprev1; m<=lprev1; m++){
+ if(agev[m][i]==0) agev[m][i]=agemax+1;
+ if(agev[m][i]==1) agev[m][i]=agemax+2;
+ 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(jk=1; jk <=nlstate ; jk++){
+ for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+ pp[jk] += freq[jk][m][i];
+ }
+ for(jk=1; jk <=nlstate ; jk++){
+ for(m=-1, pos=0; m <=0 ; m++)
+ pos += freq[jk][m][i];
+ }
+
+ 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; 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_vector(pp,1,nlstate);
+
+} /* End of Freq */
/************* Waves Concatenation ***************/
void concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm)
@@ -1492,6 +1572,12 @@ void varevsij(char fileres[], double ***
}
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+
+ if (popbased==1) {
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(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++)
@@ -1503,12 +1589,19 @@ void varevsij(char fileres[], double ***
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) {
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(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];
}
}
+
for(j=1; j<= nlstate; j++)
for(h=0; h<=nhstepm; h++){
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
@@ -1773,9 +1866,10 @@ int main()
int c, h , cpt,l;
int ju,jl, mi;
int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
- int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
-
+ int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
+ int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast;
int hstepm, nhstepm;
+
double bage, fage, age, agelim, agebase;
double ftolpl=FTOL;
double **prlim;
@@ -1852,7 +1946,25 @@ split(pathtot, path,optionfile);
fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);
fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);
-
+while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ puts(line);
+ fputs(line,ficparo);
+ }
+ ungetc(c,ficpar);
+
+ fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ puts(line);
+ fputs(line,ficparo);
+ }
+ ungetc(c,ficpar);
+
+ fscanf(ficpar,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",&fprevfore,&lprevfore,&nforecast,&mobilav);
+
covar=matrix(0,NCOVMAX,1,n);
cptcovn=0;
if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
@@ -2014,8 +2126,8 @@ split(pathtot, path,optionfile);
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[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
- }
- for (i=1; i<=imx; 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]));*/
+ }
+ for (i=1; i<=imx; 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*/
Tvar=ivector(1,15);
@@ -2206,7 +2318,7 @@ printf("Total number of individuals= %d,
/* Calculates basic frequencies. Computes observed prevalence at single age
and prints on file fileres'p'. */
- freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);
+ freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev);
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
@@ -2470,17 +2582,12 @@ ij=1;
fclose(ficgp);
chdir(path);
- free_matrix(agev,1,maxwav,1,imx);
+
free_ivector(wav,1,imx);
free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
- free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
-
- free_imatrix(s,1,maxwav+1,1,n);
-
-
+ free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
free_ivector(num,1,n);
free_vector(agedc,1,n);
- free_vector(weight,1,n);
/*free_matrix(covar,1,NCOVMAX,1,n);*/
fclose(ficparo);
fclose(ficres);
@@ -2525,7 +2632,9 @@ Interval (in months) between two waves:
- Life expectancies by age and initial health status: e%s
- Variances of life expectancies by age and initial health status: v%s
- Health expectancies with their variances: t%s
- - Standard deviation of stationary prevalences: vpl%s
",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
+ - Standard deviation of stationary prevalences: vpl%s
+ - Prevalences forecasting: f%s
+
",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
fprintf(fichtm," Graphs");
@@ -2663,7 +2772,6 @@ fclose(fichtm);
fclose(ficrespij);
- exit(0);
/*---------- Forecasting ------------------*/
strcpy(fileresf,"f");
@@ -2673,42 +2781,39 @@ fclose(fichtm);
}
printf("Computing forecasting: result on file '%s' \n", fileresf);
- /* Mobile average */
+ prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprevfore, lprevfore);
- /* for (agedeb=bage; agedeb<=fage; agedeb++)
- for (i=1; i<=nlstate;i++)
- for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++)
- printf("%f %d i=%d j1=%d\n", probs[(int)agedeb][i][cptcod],(int) agedeb,i,cptcod);*/
+ free_matrix(agev,1,maxwav,1,imx);
+ /* Mobile average */
if (cptcoveff==0) ncodemax[cptcoveff]=1;
- mobaverage= ma3x(1,130 ,1,8, 1,8);
- for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
- for (i=1; i<=nlstate;i++)
- for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
- mobaverage[(int)agedeb][i][cptcod]=0.;
-
- for (agedeb=bage+4; agedeb<=fage; agedeb++){
- for (i=1; i<=nlstate;i++){
- for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
- for (cpt=0;cpt<=4;cpt++){
- mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
+ if (mobilav==1) {
+ mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
+ for (i=1; i<=nlstate;i++)
+ for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
+ mobaverage[(int)agedeb][i][cptcod]=0.;
+
+ for (agedeb=bage+4; agedeb<=fage; agedeb++){
+ for (i=1; i<=nlstate;i++){
+ for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+ for (cpt=0;cpt<=4;cpt++){
+ mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
}
mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
+ }
}
- }
+ }
}
-/* if (cptcod==2) printf("m=%f p=%f %d age=%d ",mobaverage[(int)agedeb-2][i][cptcod],probs[(int)agedeb-cpt][i][cptcod],cpt,(int)agedeb-2);*/
-
-
stepsize=(int) (stepm+YEARM-1)/YEARM;
- if (stepm<=24) stepsize=2;
+ if (stepm<=12) stepsize=1;
agelim=AGESUP;
hstepm=stepsize*YEARM; /* Every year of age */
- hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
- hstepm=12;
+ hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */
+
k=0;
for(cptcov=1;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
@@ -2717,25 +2822,28 @@ fclose(fichtm);
for(j=1;j<=cptcoveff;j++) {
fprintf(ficresf,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
}
-
fprintf(ficresf,"******\n");
-
fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");
for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
for (agedeb=fage; agedeb>=bage; agedeb--){
fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);
+ if (mobilav==1) {
for(j=1; j<=nlstate;j++)
- fprintf(ficresf,"%.3f ",mobaverage[(int)agedeb][j][cptcod]);
+ fprintf(ficresf," %.5f ",mobaverage[(int)agedeb][j][cptcod]);
+ }
+ else {
+ for(j=1; j<=nlstate;j++)
+ fprintf(ficresf," %.5f ",probs[(int)agedeb][j][cptcod]);
+ }
+ for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");
}
- for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");
-
- for (cpt=1; cpt<=8;cpt++)
+ for (cpt=1; cpt<=nforecast;cpt++)
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
- nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
- /*printf("stepm=%d hstepm=%d nhstepm=%d \n",stepm,hstepm,nhstepm);*/
+ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
+ nhstepm = nhstepm/hstepm;
+ /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
@@ -2743,26 +2851,28 @@ fclose(fichtm);
for (h=0; h<=nhstepm; h++){
- if (h*hstepm/YEARM*stepm==cpt)
- fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);
+ if (h*hstepm/YEARM*stepm==cpt)
+ fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);
+
for(j=1; j<=nlstate+ndeath;j++) {
kk1=0.;
for(i=1; i<=nlstate;i++) {
- /* kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];*/
- kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
- }
-
- if (h*hstepm/YEARM*stepm==cpt)
- fprintf(ficresf," %.5f ", kk1);
+ if (mobilav==1)
+ kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
+ else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];
+ }
+ if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);
}
- }
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}
}
+ }
+ if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ free_imatrix(s,1,maxwav+1,1,n);
+ free_vector(weight,1,n);
fclose(ficresf);
-
/*---------- Health expectancies and variances ------------*/
strcpy(filerest,"t");
@@ -2822,6 +2932,11 @@ fclose(fichtm);
epj=vector(1,nlstate+1);
for(age=bage; age <=fage ;age++){
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+ if (popbased==1) {
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][k];
+ }
+
fprintf(ficrest," %.0f",age);
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
for(i=1, epj[j]=0.;i <=nlstate;i++) {