Moving average on cross-sectional prevalences

diff --git a/src/imach.c b/src/imach.c
index f3f7921a5eb16e4c688a0f21052f5675a90704d5..e9b172fc5967a6a952316d201c41942e432df9b6 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -73,7 +73,7 @@ int npar=NPARMAX;
 int nlstate=2; /* Number of live states */\r
 int ndeath=1; /* Number of dead states */\r
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */\r
-int popbased=0, fprev,lprev;\r
+int popbased=0;\r
 \r
 int *wav; /* Number of waves for this individuual 0 is possible */\r
 int maxwav; /* Maxim number of waves */\r
@@ -1150,7 +1150,7 @@ void lubksb(double **a, int n, int *indx, double b[])
 } \r
 \r
 /************ Frequencies ********************/\r
-void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax)\r
+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)\r
 {  /* Some frequencies */\r
  \r
   int i, m, jk, k1, i1, j1, bool, z1,z2,j;\r
@@ -1192,7 +1192,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
               bool=0;\r
         }\r
          if (bool==1) {\r
-          for(m=fprev; m<=lprev; m++){\r
+          for(m=fprev1; m<=lprev1; m++){\r
             if(agev[m][i]==0) agev[m][i]=agemax+1;\r
             if(agev[m][i]==1) agev[m][i]=agemax+2;\r
             freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];\r
@@ -1265,6 +1265,83 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
 \r
 }  /* End of Freq */\r
 \r
+/************ Prevalence ********************/\r
+void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1)\r
+{  /* Some frequencies */\r
+ \r
+  int i, m, jk, k1, i1, j1, bool, z1,z2,j;\r
+  double ***freq; /* Frequencies */\r
+  double *pp;\r
+  double pos;\r
+\r
+  pp=vector(1,nlstate);\r
+  probs= ma3x(1,130 ,1,8, 1,8);\r
+  \r
+  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);\r
+  j1=0;\r
+  \r
+  j=cptcoveff;\r
+  if (cptcovn<1) {j=1;ncodemax[1]=1;}\r
+  \r
+ for(k1=1; k1<=j;k1++){\r
+    for(i1=1; i1<=ncodemax[k1];i1++){\r
+      j1++;\r
+ \r
+      for (i=-1; i<=nlstate+ndeath; i++)  \r
+       for (jk=-1; jk<=nlstate+ndeath; jk++)  \r
+         for(m=agemin; m <= agemax+3; m++)\r
+         freq[i][jk][m]=0;\r
+      \r
+      for (i=1; i<=imx; i++) {\r
+       bool=1;\r
+       if  (cptcovn>0) {\r
+         for (z1=1; z1<=cptcoveff; z1++) \r
+           if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) \r
+             bool=0;\r
+             }\r
+       if (bool==1) {\r
+         for(m=fprev1; m<=lprev1; m++){\r
+           if(agev[m][i]==0) agev[m][i]=agemax+1;\r
+           if(agev[m][i]==1) agev[m][i]=agemax+2;\r
+           freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];\r
+           freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];\r
+         }\r
+       }\r
+      }\r
+       for(i=(int)agemin; i <= (int)agemax+3; i++){ \r
+       for(jk=1; jk <=nlstate ; jk++){\r
+         for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)\r
+           pp[jk] += freq[jk][m][i]; \r
+       }\r
+       for(jk=1; jk <=nlstate ; jk++){\r
+         for(m=-1, pos=0; m <=0 ; m++)\r
+           pos += freq[jk][m][i];\r
+       }\r
+       \r
+        for(jk=1; jk <=nlstate ; jk++){\r
+          for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)\r
+            pp[jk] += freq[jk][m][i];\r
+        }\r
+        \r
+        for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];\r
+\r
+        for(jk=1; jk <=nlstate ; jk++){           \r
+          if( i <= (int) agemax){\r
+            if(pos>=1.e-5){\r
+              probs[i][jk][j1]= pp[jk]/pos;\r
+            }\r
+          }\r
+        }\r
+        \r
+        }\r
+    }\r
+  }\r
+  \r
+  \r
+  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);\r
+  free_vector(pp,1,nlstate);\r
+  \r
+}  /* End of Freq */\r
 /************* Waves Concatenation ***************/\r
 \r
 void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)\r
@@ -1790,7 +1867,7 @@ int main()
   int ju,jl, mi;\r
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;\r
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; \r
-  int mobilav=0, fprevfore=1, lprevfore=1;\r
+  int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast;\r
   int hstepm, nhstepm;\r
 \r
   double bage, fage, age, agelim, agebase;\r
@@ -1886,7 +1963,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   }\r
   ungetc(c,ficpar);\r
   \r
-  fscanf(ficpar,"fprevalence=%d lprevalence=%d mob_average=%d\n",&fprevfore,&lprevfore,&mobilav);\r
+  fscanf(ficpar,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",&fprevfore,&lprevfore,&nforecast,&mobilav);\r
  \r
   covar=matrix(0,NCOVMAX,1,n); \r
   cptcovn=0; \r
@@ -2241,7 +2318,7 @@ printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx,
     \r
    /* Calculates basic frequencies. Computes observed prevalence at single age\r
        and prints on file fileres'p'. */\r
-  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);\r
+  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev);\r
 \r
     pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */\r
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */\r
@@ -2505,17 +2582,12 @@ ij=1;
   fclose(ficgp);\r
    \r
 chdir(path);\r
-    free_matrix(agev,1,maxwav,1,imx);\r
+   \r
     free_ivector(wav,1,imx);\r
     free_imatrix(dh,1,lastpass-firstpass+1,1,imx);\r
-    free_imatrix(mw,1,lastpass-firstpass+1,1,imx);\r
-    \r
-    free_imatrix(s,1,maxwav+1,1,n);\r
-    \r
-    \r
+    free_imatrix(mw,1,lastpass-firstpass+1,1,imx);   \r
     free_ivector(num,1,n);\r
     free_vector(agedc,1,n);\r
-    free_vector(weight,1,n);\r
     /*free_matrix(covar,1,NCOVMAX,1,n);*/\r
     fclose(ficparo);\r
     fclose(ficres);\r
@@ -2709,12 +2781,10 @@ fclose(fichtm);
   }\r
   printf("Computing forecasting: result on file '%s' \n", fileresf);\r
 \r
-  /* Mobile average */\r
+  prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprevfore, lprevfore);\r
 \r
-  /* for (agedeb=bage; agedeb<=fage; agedeb++)\r
-    for (i=1; i<=nlstate;i++)\r
-      for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++)\r
-      printf("%f %d i=%d j1=%d\n", probs[(int)agedeb][i][cptcod],(int) agedeb,i,cptcod);*/\r
+ free_matrix(agev,1,maxwav,1,imx);\r
+  /* Mobile average */\r
 \r
   if (cptcoveff==0) ncodemax[cptcoveff]=1;\r
 \r
@@ -2738,12 +2808,12 @@ fclose(fichtm);
   }\r
 \r
   stepsize=(int) (stepm+YEARM-1)/YEARM;\r
-  if (stepm<=24) stepsize=2;\r
+  if (stepm<=12) stepsize=1;\r
 \r
   agelim=AGESUP;\r
   hstepm=stepsize*YEARM; /* Every year of age */\r
-  hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ \r
-  hstepm=12;\r
+  hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */ \r
+\r
    k=0;\r
   for(cptcov=1;cptcov<=i1;cptcov++){\r
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
@@ -2752,9 +2822,7 @@ fclose(fichtm);
       for(j=1;j<=cptcoveff;j++) {\r
        fprintf(ficresf,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
       }\r
-      \r
       fprintf(ficresf,"******\n");\r
-\r
       fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");\r
       for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);\r
 \r
@@ -2762,21 +2830,20 @@ fclose(fichtm);
        fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);\r
        if (mobilav==1) {\r
        for(j=1; j<=nlstate;j++) \r
-         fprintf(ficresf,"%.5f ",mobaverage[(int)agedeb][j][cptcod]);\r
+         fprintf(ficresf," %.5f ",mobaverage[(int)agedeb][j][cptcod]);\r
        }\r
        else {\r
          for(j=1; j<=nlstate;j++) \r
-         fprintf(ficresf,"%.5f ",probs[(int)agedeb][j][cptcod]);\r
-       }\r
-       \r
-      for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");\r
+         fprintf(ficresf," %.5f ",probs[(int)agedeb][j][cptcod]);\r
+       }  \r
+      for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");\r
       }\r
-      for (cpt=1; cpt<=NCOVMAX;cpt++)   \r
+      for (cpt=1; cpt<=nforecast;cpt++)   \r
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */\r
        \r
-       nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ \r
-       nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */\r
-       /*printf("stepm=%d hstepm=%d nhstepm=%d \n",stepm,hstepm,nhstepm);*/\r
+       nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
+       nhstepm = nhstepm/hstepm; \r
+       /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/\r
 \r
        p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
        oldm=oldms;savm=savms;\r
@@ -2784,8 +2851,9 @@ fclose(fichtm);
                \r
        for (h=0; h<=nhstepm; h++){\r
        \r
-         if (h*hstepm/YEARM*stepm==cpt)\r
- fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);\r
+        if (h*hstepm/YEARM*stepm==cpt)\r
+           fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);\r
+        \r
          \r
          for(j=1; j<=nlstate+ndeath;j++) {\r
            kk1=0.;\r
@@ -2794,7 +2862,7 @@ fclose(fichtm);
              kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];\r
              else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];\r
            }    \r
-           if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);\r
+         if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);\r
          }\r
        }\r
        free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
@@ -2802,6 +2870,8 @@ fclose(fichtm);
     }\r
   }\r
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+  free_imatrix(s,1,maxwav+1,1,n);\r
+  free_vector(weight,1,n);\r
   fclose(ficresf);\r
   /*---------- Health expectancies and variances ------------*/\r
 \r