--- imach/src/imach.c 2000/12/28 18:49:56 1.1 +++ imach/src/imach.c 2009/07/06 08:22:05 1.132 @@ -1,37 +1,308 @@ - -/*********************** Imach ************************************** - This program computes Healthy Life Expectancies from cross-longitudinal - data. Cross-longitudinal consist in a first survey ("cross") where - individuals from different ages are interviewed on their health status - or degree of disability. At least a second wave of interviews - ("longitudinal") should measure each new individual health status. - Health expectancies are computed from the transistions observed between - waves and are computed for each degree of severity of disability (number - of life states). More degrees you consider, more time is necessary to - reach the Maximum Likekilhood of the parameters involved in the model. - The simplest model is the multinomial logistic model where pij is - the probabibility to be observed in state j at the second wave 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 'age' is age and 'sex' - is a covariate. If you want to have a more complex model than "constant and - age", you should modify the program where the markup - *Covariates have to be included here again* invites you to do it. - More covariates you add, less is the speed of the convergence. - - The advantage that this computer programme claims, comes from that if the - delay between waves is not identical for each individual, or if some - individual missed an interview, the information is not rounded or lost, but - taken into account using an interpolation or extrapolation. - hPijx is the probability to be - observed in state i at age x+h conditional to the observed state i at age - x. The delay 'h' can be split into an exact number (nh*stepm) of - unobserved intermediate states. This elementary transition (by month or - quarter trimester, semester or year) is model as a multinomial logistic. - The hPx matrix is simply the matrix product of nh*stepm elementary matrices - and the contribution of each individual to the likelihood is simply hPijx. +/* $Id: imach.c,v 1.132 2009/07/06 08:22:05 brouard Exp $ + $State: Exp $ + $Log: imach.c,v $ + Revision 1.132 2009/07/06 08:22:05 brouard + Many tings + + Revision 1.131 2009/06/20 16:22:47 brouard + Some dimensions resccaled + + Revision 1.130 2009/05/26 06:44:34 brouard + (Module): Max Covariate is now set to 20 instead of 8. A + lot of cleaning with variables initialized to 0. Trying to make + V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better. + + Revision 1.129 2007/08/31 13:49:27 lievre + Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting + + Revision 1.128 2006/06/30 13:02:05 brouard + (Module): Clarifications on computing e.j + + Revision 1.127 2006/04/28 18:11:50 brouard + (Module): Yes the sum of survivors was wrong since + imach-114 because nhstepm was no more computed in the age + loop. Now we define nhstepma in the age loop. + (Module): In order to speed up (in case of numerous covariates) we + compute health expectancies (without variances) in a first step + and then all the health expectancies with variances or standard + deviation (needs data from the Hessian matrices) which slows the + computation. + In the future we should be able to stop the program is only health + expectancies and graph are needed without standard deviations. + + Revision 1.126 2006/04/28 17:23:28 brouard + (Module): Yes the sum of survivors was wrong since + imach-114 because nhstepm was no more computed in the age + loop. Now we define nhstepma in the age loop. + Version 0.98h + + Revision 1.125 2006/04/04 15:20:31 lievre + Errors in calculation of health expectancies. Age was not initialized. + Forecasting file added. + + Revision 1.124 2006/03/22 17:13:53 lievre + Parameters are printed with %lf instead of %f (more numbers after the comma). + The log-likelihood is printed in the log file + + Revision 1.123 2006/03/20 10:52:43 brouard + * imach.c (Module):
=(p+1))(v[j-p-1] = t[j]);
}
}
@@ -166,7 +654,7 @@ void nrerror(char error_text[])
{
fprintf(stderr,"ERREUR ...\n");
fprintf(stderr,"%s\n",error_text);
- exit(1);
+ exit(EXIT_FAILURE);
}
/*********************** vector *******************/
double *vector(int nl, int nh)
@@ -198,6 +686,21 @@ void free_ivector(int *v, long nl, long
free((FREE_ARG)(v+nl-NR_END));
}
+/************************lvector *******************************/
+long *lvector(long nl,long nh)
+{
+ long *v;
+ v=(long *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(long)));
+ if (!v) nrerror("allocation failure in ivector");
+ return v-nl+NR_END;
+}
+
+/******************free lvector **************************/
+void free_lvector(long *v, long nl, long nh)
+{
+ free((FREE_ARG)(v+nl-NR_END));
+}
+
/******************* imatrix *******************************/
int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
@@ -252,6 +755,8 @@ double **matrix(long nrl, long nrh, long
for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
return m;
+ /* print *(*(m+1)+70) or print m[1][70]; print m+1 or print &(m[1])
+ */
}
/*************************free matrix ************************/
@@ -291,7 +796,10 @@ double ***ma3x(long nrl, long nrh, long
for (j=ncl+1; j<=nch; j++)
m[i][j]=m[i][j-1]+nlay;
}
- return m;
+ return m;
+ /* gdb: p *(m+1) <=> p m[1] and p (m+1) <=> p (m+1) <=> p &(m[1])
+ &(m[i][j][k]) <=> *((*(m+i) + j)+k)
+ */
}
/*************************free ma3x ************************/
@@ -302,6 +810,41 @@ void free_ma3x(double ***m, long nrl, lo
free((FREE_ARG)(m+nrl-NR_END));
}
+/*************** function subdirf ***********/
+char *subdirf(char fileres[])
+{
+ /* Caution optionfilefiname is hidden */
+ strcpy(tmpout,optionfilefiname);
+ strcat(tmpout,"/"); /* Add to the right */
+ strcat(tmpout,fileres);
+ return tmpout;
+}
+
+/*************** function subdirf2 ***********/
+char *subdirf2(char fileres[], char *preop)
+{
+
+ /* Caution optionfilefiname is hidden */
+ strcpy(tmpout,optionfilefiname);
+ strcat(tmpout,"/");
+ strcat(tmpout,preop);
+ strcat(tmpout,fileres);
+ return tmpout;
+}
+
+/*************** function subdirf3 ***********/
+char *subdirf3(char fileres[], char *preop, char *preop2)
+{
+
+ /* Caution optionfilefiname is hidden */
+ strcpy(tmpout,optionfilefiname);
+ strcat(tmpout,"/");
+ strcat(tmpout,preop);
+ strcat(tmpout,preop2);
+ strcat(tmpout,fileres);
+ return tmpout;
+}
+
/***************** f1dim *************************/
extern int ncom;
extern double *pcom,*xicom;
@@ -339,8 +882,10 @@ double brent(double ax, double bx, doubl
tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/
printf(".");fflush(stdout);
+ fprintf(ficlog,".");fflush(ficlog);
#ifdef DEBUG
printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
+ fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
/* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */
#endif
if (fabs(x-xm) <= (tol2-0.5*(b-a))){
@@ -465,6 +1010,7 @@ void linmin(double p[], double xi[], int
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
#ifdef DEBUG
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
+ fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
#endif
for (j=1;j<=n;j++) {
xi[j] *= xmin;
@@ -474,19 +1020,31 @@ void linmin(double p[], double xi[], int
free_vector(pcom,1,n);
}
+char *asc_diff_time(long time_sec, char ascdiff[])
+{
+ long sec_left, days, hours, minutes;
+ days = (time_sec) / (60*60*24);
+ sec_left = (time_sec) % (60*60*24);
+ hours = (sec_left) / (60*60) ;
+ sec_left = (sec_left) %(60*60);
+ minutes = (sec_left) /60;
+ sec_left = (sec_left) % (60);
+ sprintf(ascdiff,"%d day(s) %d hour(s) %d minute(s) %d second(s)",days, hours, minutes, sec_left);
+ return ascdiff;
+}
+
/*************** powell ************************/
void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,
double (*func)(double []))
-
{
-
-
void linmin(double p[], double xi[], int n, double *fret,
double (*func)(double []));
int i,ibig,j;
double del,t,*pt,*ptt,*xit;
double fp,fptt;
double *xits;
+ int niterf, itmp;
+
pt=vector(1,n);
ptt=vector(1,n);
xit=vector(1,n);
@@ -497,17 +1055,50 @@ void powell(double p[], double **xi, int
fp=(*fret);
ibig=0;
del=0.0;
- printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
- for (i=1;i<=n;i++)
+ last_time=curr_time;
+ (void) gettimeofday(&curr_time,&tzp);
+ printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout);
+ fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fflush(ficlog);
+/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */
+ for (i=1;i<=n;i++) {
printf(" %d %.12f",i, p[i]);
+ fprintf(ficlog," %d %.12lf",i, p[i]);
+ fprintf(ficrespow," %.12lf", p[i]);
+ }
printf("\n");
+ fprintf(ficlog,"\n");
+ fprintf(ficrespow,"\n");fflush(ficrespow);
+ if(*iter <=3){
+ tm = *localtime(&curr_time.tv_sec);
+ strcpy(strcurr,asctime(&tm));
+/* asctime_r(&tm,strcurr); */
+ forecast_time=curr_time;
+ itmp = strlen(strcurr);
+ if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */
+ strcurr[itmp-1]='\0';
+ printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
+ fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
+ for(niterf=10;niterf<=30;niterf+=10){
+ forecast_time.tv_sec=curr_time.tv_sec+(niterf-*iter)*(curr_time.tv_sec-last_time.tv_sec);
+ tmf = *localtime(&forecast_time.tv_sec);
+/* asctime_r(&tmf,strfor); */
+ strcpy(strfor,asctime(&tmf));
+ itmp = strlen(strfor);
+ if(strfor[itmp-1]=='\n')
+ strfor[itmp-1]='\0';
+ printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+ fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+ }
+ }
for (i=1;i<=n;i++) {
for (j=1;j<=n;j++) xit[j]=xi[j][i];
fptt=(*fret);
#ifdef DEBUG
printf("fret=%lf \n",*fret);
+ fprintf(ficlog,"fret=%lf \n",*fret);
#endif
printf("%d",i);fflush(stdout);
+ fprintf(ficlog,"%d",i);fflush(ficlog);
linmin(p,xit,n,fret,func);
if (fabs(fptt-(*fret)) > del) {
del=fabs(fptt-(*fret));
@@ -515,13 +1106,18 @@ void powell(double p[], double **xi, int
}
#ifdef DEBUG
printf("%d %.12e",i,(*fret));
+ fprintf(ficlog,"%d %.12e",i,(*fret));
for (j=1;j<=n;j++) {
xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
printf(" x(%d)=%.12e",j,xit[j]);
+ fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
}
- for(j=1;j<=n;j++)
+ for(j=1;j<=n;j++) {
printf(" p=%.12e",p[j]);
+ fprintf(ficlog," p=%.12e",p[j]);
+ }
printf("\n");
+ fprintf(ficlog,"\n");
#endif
}
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
@@ -530,15 +1126,21 @@ void powell(double p[], double **xi, int
k[0]=1;
k[1]=-1;
printf("Max: %.12e",(*func)(p));
- for (j=1;j<=n;j++)
+ fprintf(ficlog,"Max: %.12e",(*func)(p));
+ for (j=1;j<=n;j++) {
printf(" %.12e",p[j]);
+ fprintf(ficlog," %.12e",p[j]);
+ }
printf("\n");
+ fprintf(ficlog,"\n");
for(l=0;l<=1;l++) {
for (j=1;j<=n;j++) {
ptt[j]=p[j]+(p[j]-pt[j])*k[l];
printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
+ fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
}
printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
+ fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
}
#endif
@@ -566,18 +1168,22 @@ void powell(double p[], double **xi, int
}
#ifdef DEBUG
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
- for(j=1;j<=n;j++)
+ fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+ for(j=1;j<=n;j++){
printf(" %.12e",xit[j]);
+ fprintf(ficlog," %.12e",xit[j]);
+ }
printf("\n");
+ fprintf(ficlog,"\n");
#endif
- }
+ }
}
}
}
-/**** Prevalence limit ****************/
+/**** Prevalence limit (stable or period prevalence) ****************/
-double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl)
+double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
{
/* Computes the prevalence limit in each live state at age x by left multiplying the unit
matrix by transitions matrix until convergence is reached */
@@ -585,7 +1191,7 @@ double **prevalim(double **prlim, int nl
int i, ii,j,k;
double min, max, maxmin, maxmax,sumnew=0.;
double **matprod2();
- double **out, cov[NCOVMAX], **pmij();
+ double **out, cov[NCOVMAX+1], **pmij();
double **newm;
double agefin, delaymax=50 ; /* Max number of years to converge */
@@ -593,15 +1199,28 @@ double **prevalim(double **prlim, int nl
for (j=1;j<=nlstate+ndeath;j++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);
}
- /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+
+ cov[1]=1.;
+
+ /* Even if hstepm = 1, at least one multiplication by the unit matrix */
for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
newm=savm;
/* Covariates have to be included here again */
- cov[1]=1.;
- cov[2]=agefin;
- out=matprod2(newm, pmij(pmmij,cov,ncov,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
-/* printf("age=%f agefin=%f po=%f pn=%f\n",age,agefin,oldm[1][1],newm[1][1]);*/
-
+ cov[2]=agefin;
+
+ for (k=1; k<=cptcovn;k++) {
+ cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[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++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
+ for (k=1; k<=cptcovprod;k++)
+ cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+
+ /*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 \n",ij, cov[3]);*/
+ out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+
savm=oldm;
oldm=newm;
maxmax=0.;
@@ -624,63 +1243,70 @@ double **prevalim(double **prlim, int nl
}
}
-/*************** transition probabilities **********/
+/*************** transition probabilities ***************/
-double **pmij(double **ps, double *cov, int ncov, double *x, int nlstate )
+double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
{
double s1, s2;
/*double t34;*/
int i,j,j1, nc, ii, jj;
for(i=1; i<= nlstate; i++){
- for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
- }
- ps[i][j]=s2;
- }
- }
- for(i=1; i<= nlstate; i++){
- s1=0;
- for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */
+ }
+ ps[i][j]=s2;
+ }
}
- }
-
- /* for(ii=1; ii<= nlstate+ndeath; ii++){
- for(jj=1; jj<= nlstate+ndeath; jj++){
- printf("%lf ",ps[ii][jj]);
- }
- printf("\n ");
+ /*ps[3][2]=1;*/
+
+ for(i=1; i<= nlstate; i++){
+ s1=0;
+ for(j=1; j 1 the results are less biased than in previous versions.
+ */
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ /* bias bh is positive if real duration
+ * is higher than the multiple of stepm and negative otherwise.
+ */
+ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
+ if( s2 > nlstate){
+ /* i.e. if s2 is a death state and if the date of death is known
+ then the contribution to the likelihood is the probability to
+ die between last step unit time and current step unit time,
+ which is also equal to probability to die before dh
+ minus probability to die before dh-stepm .
+ In version up to 0.92 likelihood was computed
+ as if date of death was unknown. Death was treated as any other
+ health state: the date of the interview describes the actual state
+ and not the date of a change in health state. The former idea was
+ to consider that at each interview the state was recorded
+ (healthy, disable or death) and IMaCh was corrected; but when we
+ introduced the exact date of death then we should have modified
+ the contribution of an exact death to the likelihood. This new
+ contribution is smaller and very dependent of the step unit
+ stepm. It is no more the probability to die between last interview
+ and month of death but the probability to survive from last
+ interview up to one month before death multiplied by the
+ probability to die within a month. Thanks to Chris
+ Jackson for correcting this bug. Former versions increased
+ mortality artificially. The bad side is that we add another loop
+ which slows down the processing. The difference can be up to 10%
+ lower mortality.
+ */
+ lli=log(out[s1][s2] - savm[s1][s2]);
+
+
+ } else if (s2==-2) {
+ for (j=1,survp=0. ; j<=nlstate; j++)
+ survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ /*survp += out[s1][j]; */
+ lli= log(survp);
+ }
+
+ else if (s2==-4) {
+ for (j=3,survp=0. ; j<=nlstate; j++)
+ survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ lli= log(survp);
+ }
+
+ else if (s2==-5) {
+ for (j=1,survp=0. ; j<=2; j++)
+ survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ lli= log(survp);
+ }
+
+ else{
+ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+ /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
+ }
+ /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
+ /*if(lli ==000.0)*/
+ /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ } /* end of wave */
+ } /* end of individual */
+ } else if(mle==2){
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<=dh[mi][i]; d++){
+ newm=savm;
+ cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
+
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ } /* end of wave */
+ } /* end of individual */
+ } else if(mle==3){ /* exponential inter-extrapolation */
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d ");
+
+ m=cptcoveff;
+ if (cptcovn < 1) {m=1;ncodemax[1]=1;}
+
+ jj1=0;
+ for(k1=1; k1<=m;k1++){
+ for(i1=1; i1<=ncodemax[k1];i1++){
+ jj1++;
+ if (cptcovn > 0) {
+ fprintf(fichtm," ");
+
+ m=cptcoveff;
+ if (cptcovn < 1) {m=1;ncodemax[1]=1;}
+
+ jj1=0;
+ for(k1=1; k1<=m;k1++){
+ for(i1=1; i1<=ncodemax[k1];i1++){
+ jj1++;
+ if (cptcovn > 0) {
+ fprintf(fichtm,"
File of contributions to the likelihood: %s
\n",subdirf(fileresilk),subdirf(fileresilk));
+ fflush(fichtm);
+ }
+ return;
+}
+
+
/*********** Maximum Likelihood Estimation ***************/
-void mlikeli(FILE *ficres,double p[], int npar, int ncov, int nlstate, double ftol, double (*func)(double []))
+void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))
{
int i,j, iter;
- double **xi,*delti;
+ double **xi;
double fret;
+ double fretone; /* Only one call to likelihood */
+ /* char filerespow[FILENAMELENGTH];*/
xi=matrix(1,npar,1,npar);
for (i=1;i<=npar;i++)
for (j=1;j<=npar;j++)
xi[i][j]=(i==j ? 1.0 : 0.0);
- printf("Powell\n");
+ printf("Powell\n"); fprintf(ficlog,"Powell\n");
+ strcpy(filerespow,"pow");
+ strcat(filerespow,fileres);
+ if((ficrespow=fopen(filerespow,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", filerespow);
+ fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
+ }
+ fprintf(ficrespow,"# Powell\n# iter -2*LL");
+ for (i=1;i<=nlstate;i++)
+ for(j=1;j<=nlstate+ndeath;j++)
+ if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
+ fprintf(ficrespow,"\n");
+
powell(p,xi,npar,ftol,&iter,&fret,func);
- printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
- fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f ",iter,func(p));
+ free_matrix(xi,1,npar,1,npar);
+ fclose(ficrespow);
+ printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
+ fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
+ fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
}
@@ -826,33 +1813,42 @@ void hesscov(double **matcov, double p[]
int i, j,jk;
int *indx;
- double hessii(double p[], double delta, int theta, double delti[]);
- double hessij(double p[], double delti[], int i, int j);
+ double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar);
+ double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar);
void lubksb(double **a, int npar, int *indx, double b[]) ;
void ludcmp(double **a, int npar, int *indx, double *d) ;
-
-
+ double gompertz(double p[]);
hess=matrix(1,npar,1,npar);
printf("\nCalculation of the hessian matrix. Wait...\n");
+ fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n");
for (i=1;i<=npar;i++){
printf("%d",i);fflush(stdout);
- hess[i][i]=hessii(p,ftolhess,i,delti);
- /*printf(" %f ",p[i]);*/
+ fprintf(ficlog,"%d",i);fflush(ficlog);
+
+ hess[i][i]=hessii(p,ftolhess,i,delti,func,npar);
+
+ /* printf(" %f ",p[i]);
+ printf(" %lf %lf %lf",hess[i][i],ftolhess,delti[i]);*/
}
-
+
for (i=1;i<=npar;i++) {
for (j=1;j<=npar;j++) {
if (j>i) {
printf(".%d%d",i,j);fflush(stdout);
- hess[i][j]=hessij(p,delti,i,j);
- hess[j][i]=hess[i][j];
+ fprintf(ficlog,".%d%d",i,j);fflush(ficlog);
+ hess[i][j]=hessij(p,delti,i,j,func,npar);
+
+ hess[j][i]=hess[i][j];
+ /*printf(" %lf ",hess[i][j]);*/
}
}
}
printf("\n");
+ fprintf(ficlog,"\n");
printf("\nInverting the hessian to get the covariance matrix. Wait...\n");
+ fprintf(ficlog,"\nInverting the hessian to get the covariance matrix. Wait...\n");
a=matrix(1,npar,1,npar);
y=matrix(1,npar,1,npar);
@@ -872,11 +1868,14 @@ void hesscov(double **matcov, double p[]
}
printf("\n#Hessian matrix#\n");
+ fprintf(ficlog,"\n#Hessian matrix#\n");
for (i=1;i<=npar;i++) {
for (j=1;j<=npar;j++) {
printf("%.3e ",hess[i][j]);
+ fprintf(ficlog,"%.3e ",hess[i][j]);
}
printf("\n");
+ fprintf(ficlog,"\n");
}
/* Recompute Inverse */
@@ -893,8 +1892,10 @@ void hesscov(double **matcov, double p[]
for (i=1;i<=npar;i++){
y[i][j]=x[i];
printf("%.3e ",y[i][j]);
+ fprintf(ficlog,"%.3e ",y[i][j]);
}
printf("\n");
+ fprintf(ficlog,"\n");
}
*/
@@ -908,14 +1909,14 @@ void hesscov(double **matcov, double p[]
}
/*************** hessian matrix ****************/
-double hessii( double x[], double delta, int theta, double delti[])
+double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar)
{
int i;
int l=1, lmax=20;
double k1,k2;
- double p2[NPARMAX+1];
+ double p2[MAXPARM+1]; /* identical to x */
double res;
- double delt, delts, nkhi=10.,nkhif=1., khi=1.e-4;
+ double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4;
double fx;
int k=0,kmax=10;
double l1;
@@ -934,8 +1935,9 @@ double hessii( double x[], double delta,
/*res= (k1-2.0*fx+k2)/delt/delt; */
res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */
-#ifdef DEBUG
+#ifdef DEBUGHESS
printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);
+ fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);
#endif
/*if(fabs(k1-2.0*fx+k2) <1.e-13){ */
if((k1 Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)
%s
\n",digitp);
+/* } */
+ varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+ pstamp(ficresvij);
+ fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are ");
+ if(popbased==1)
+ fprintf(ficresvij,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d\n",mobilav);
+ else
+ fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");
fprintf(ficresvij,"# Age");
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++)
- fprintf(ficresvij," Cov(e%1d, e%1d)",i,j);
+ fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
fprintf(ficresvij,"\n");
xp=vector(1,npar);
dnewm=matrix(1,nlstate,1,npar);
doldm=matrix(1,nlstate,1,nlstate);
-
- hstepm=1*YEARM; /* Every year of age */
- hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */
+ dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);
+ doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+
+ gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);
+ gpp=vector(nlstate+1,nlstate+ndeath);
+ gmp=vector(nlstate+1,nlstate+ndeath);
+ trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
+
+ if(estepm < stepm){
+ printf ("Problem %d lower than %d\n",estepm, stepm);
+ }
+ else hstepm=estepm;
+ /* For example we decided to compute the life expectancy with the smallest unit */
+ /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
+ nhstepm is the number of hstepm from age to agelim
+ nstepm is the number of stepm from age to agelin.
+ Look at function hpijx to understand why (it is linked to memory size questions) */
+ /* We decided (b) to get a life expectancy respecting the most precise curvature of the
+ survival function given by stepm (the optimization length). Unfortunately it
+ means that if the survival funtion is printed every two years of age and if
+ you sum them up and add 1 year (area under the trapezoids) you won't get the same
+ results. So we changed our mind and took the option of the best precision.
+ */
+ hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
agelim = AGESUP;
for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
- nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- if (stepm >= YEARM) hstepm=1;
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+ nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
gp=matrix(0,nhstepm,1,nlstate);
gm=matrix(0,nhstepm,1,nlstate);
+
for(theta=1; theta <=npar; theta++){
- for(i=1; i<=npar; i++){ /* Computes gradient */
+ for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+
+ if (popbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][ij];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][ij];
+ }
+ }
+
for(j=1; j<= nlstate; j++){
for(h=0; h<=nhstepm; h++){
for(i=1, gp[h][j]=0.;i<=nlstate;i++)
gp[h][j] += prlim[i][i]*p3mat[i][j][h];
}
}
-
- for(i=1; i<=npar; i++) /* Computes gradient */
+ /* This for computing probability of death (h=1 means
+ computed over hstepm matrices product = hstepm*stepm months)
+ as a weighted average of prlim.
+ */
+ for(j=nlstate+1;j<=nlstate+ndeath;j++){
+ for(i=1,gpp[j]=0.; i<= nlstate; i++)
+ gpp[j] += prlim[i][i]*p3mat[i][j][1];
+ }
+ /* end probability of death */
+
+ for(i=1; i<=npar; i++) /* Computes gradient x - delta */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);
- for(j=1; j<= nlstate; j++){
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+
+ if (popbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][ij];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][ij];
+ }
+ }
+
+ for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */
for(h=0; h<=nhstepm; h++){
for(i=1, gm[h][j]=0.;i<=nlstate;i++)
gm[h][j] += prlim[i][i]*p3mat[i][j][h];
}
}
- for(j=1; j<= nlstate; j++)
+ /* This for computing probability of death (h=1 means
+ computed over hstepm matrices product = hstepm*stepm months)
+ as a weighted average of prlim.
+ */
+ for(j=nlstate+1;j<=nlstate+ndeath;j++){
+ for(i=1,gmp[j]=0.; i<= nlstate; i++)
+ gmp[j] += prlim[i][i]*p3mat[i][j][1];
+ }
+ /* end probability of death */
+
+ for(j=1; j<= nlstate; j++) /* vareij */
for(h=0; h<=nhstepm; h++){
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
}
+
+ for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
+ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
+ }
+
} /* End theta */
- trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);
+ trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
- for(h=0; h<=nhstepm; h++)
+ for(h=0; h<=nhstepm; h++) /* veij */
for(j=1; j<=nlstate;j++)
for(theta=1; theta <=npar; theta++)
trgradg[h][j][theta]=gradg[h][theta][j];
+ for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
+ for(theta=1; theta <=npar; theta++)
+ trgradgp[j][theta]=gradgp[theta][j];
+
+
+ hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
for(i=1;i<=nlstate;i++)
for(j=1;j<=nlstate;j++)
vareij[i][j][(int)age] =0.;
+
for(h=0;h<=nhstepm;h++){
for(k=0;k<=nhstepm;k++){
matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
for(i=1;i<=nlstate;i++)
for(j=1;j<=nlstate;j++)
- vareij[i][j][(int)age] += doldm[i][j];
+ vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
+ }
+ }
+
+ /* pptj */
+ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
+ matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
+ for(j=nlstate+1;j<=nlstate+ndeath;j++)
+ for(i=nlstate+1;i<=nlstate+ndeath;i++)
+ varppt[j][i]=doldmp[j][i];
+ /* end ppptj */
+ /* x centered again */
+ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
+ prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);
+
+ if (popbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][ij];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][ij];
}
}
- h=1;
- if (stepm >= YEARM) h=stepm/YEARM;
+
+ /* This for computing probability of death (h=1 means
+ computed over hstepm (estepm) matrices product = hstepm*stepm months)
+ as a weighted average of prlim.
+ */
+ for(j=nlstate+1;j<=nlstate+ndeath;j++){
+ for(i=1,gmp[j]=0.;i<= nlstate; i++)
+ gmp[j] += prlim[i][i]*p3mat[i][j][1];
+ }
+ /* end probability of death */
+
+ fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
+ for(j=nlstate+1; j<=(nlstate+ndeath);j++){
+ fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));
+ for(i=1; i<=nlstate;i++){
+ fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
+ }
+ }
+ fprintf(ficresprobmorprev,"\n");
+
fprintf(ficresvij,"%.0f ",age );
for(i=1; i<=nlstate;i++)
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");
free_matrix(gp,0,nhstepm,1,nlstate);
@@ -1392,40 +3112,55 @@ void varevsij(char fileres[], double ***
free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
} /* End age */
- fclose(ficresvij);
- free_vector(xp,1,npar);
- free_matrix(doldm,1,nlstate,1,npar);
- free_matrix(dnewm,1,nlstate,1,nlstate);
+ free_vector(gpp,nlstate+1,nlstate+ndeath);
+ free_vector(gmp,nlstate+1,nlstate+ndeath);
+ free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
+ free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
+ fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65");
+ /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
+ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
+/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
+/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
+ fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l 1 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l 2 ",subdirf(fileresprobmorprev));
+ fprintf(fichtm,"\n
File (multiple files are possible if covariates are present): %s\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
+ fprintf(fichtm,"\n
Probability is computed over estepm=%d months.
\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
+ /* fprintf(fichtm,"\n
Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year
\n", stepm,YEARM,digitp,digit);
+*/
+/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit); */
+ fprintf(ficgp,"\nset out \"%s%s.png\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
-}
+ free_vector(xp,1,npar);
+ free_matrix(doldm,1,nlstate,1,nlstate);
+ free_matrix(dnewm,1,nlstate,1,npar);
+ free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+ free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
+ free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+ if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ fclose(ficresprobmorprev);
+ fflush(ficgp);
+ fflush(fichtm);
+} /* end varevsij */
/************ Variance of prevlim ******************/
-void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl)
+void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[])
{
- /* Variance of health expectancies */
- /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
+ /* Variance of prevalence limit */
+ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
double **newm;
double **dnewm,**doldm;
int i, j, nhstepm, hstepm;
- int k;
- FILE *ficresvpl;
- char fileresvpl[FILENAMELENGTH];
+ int k, cptcode;
double *xp;
double *gp, *gm;
double **gradg, **trgradg;
double age,agelim;
int theta;
-
- strcpy(fileresvpl,"vpl");
- strcat(fileresvpl,fileres);
- if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
- printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
- exit(0);
- }
- printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
-
-
- fprintf(ficresvpl,"# Standard deviation of prevalences limit\n");
+
+ pstamp(ficresvpl);
+ fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
fprintf(ficresvpl,"# Age");
for(i=1; i<=nlstate;i++)
fprintf(ficresvpl," %1d-%1d",i,i);
@@ -1450,13 +3185,13 @@ void varprevlim(char fileres[], double *
for(i=1; i<=npar; i++){ /* Computes gradient */
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
for(i=1;i<=nlstate;i++)
gp[i] = prlim[i][i];
for(i=1; i<=npar; i++) /* Computes gradient */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
for(i=1;i<=nlstate;i++)
gm[i] = prlim[i][i];
@@ -1486,435 +3221,2403 @@ void varprevlim(char fileres[], double *
free_matrix(gradg,1,npar,1,nlstate);
free_matrix(trgradg,1,nlstate,1,npar);
} /* End age */
- fclose(ficresvpl);
+
free_vector(xp,1,npar);
free_matrix(doldm,1,nlstate,1,npar);
free_matrix(dnewm,1,nlstate,1,nlstate);
}
-
-
-/***********************************************/
-/**************** Main Program *****************/
-/***********************************************/
-
-/*int main(int argc, char *argv[])*/
-int main()
+/************ Variance of one-step probabilities ******************/
+void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[])
{
+ int i, j=0, i1, k1, l1, t, tj;
+ int k2, l2, j1, z1;
+ int k=0,l, cptcode;
+ int first=1, first1;
+ double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
+ double **dnewm,**doldm;
+ double *xp;
+ double *gp, *gm;
+ double **gradg, **trgradg;
+ double **mu;
+ double age,agelim, cov[NCOVMAX];
+ double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
+ int theta;
+ char fileresprob[FILENAMELENGTH];
+ char fileresprobcov[FILENAMELENGTH];
+ char fileresprobcor[FILENAMELENGTH];
+
+ double ***varpij;
+
+ strcpy(fileresprob,"prob");
+ strcat(fileresprob,fileres);
+ if((ficresprob=fopen(fileresprob,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", fileresprob);
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
+ }
+ strcpy(fileresprobcov,"probcov");
+ strcat(fileresprobcov,fileres);
+ if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", fileresprobcov);
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
+ }
+ strcpy(fileresprobcor,"probcor");
+ strcat(fileresprobcor,fileres);
+ if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", fileresprobcor);
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
+ }
+ printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
+ fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
+ printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
+ fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
+ printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
+ fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
+ pstamp(ficresprob);
+ fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
+ fprintf(ficresprob,"# Age");
+ pstamp(ficresprobcov);
+ fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
+ fprintf(ficresprobcov,"# Age");
+ pstamp(ficresprobcor);
+ fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
+ fprintf(ficresprobcor,"# Age");
- int i,j, k, n=MAXN,iter,m,size;
- double agedeb, agefin,hf;
- double agemin=1.e20, agemax=-1.e20;
-
- double fret;
- double **xi,tmp,delta;
- double dum; /* Dummy variable */
- double ***p3mat;
- int *indx;
- char line[MAXLINE], linepar[MAXLINE];
- char title[MAXLINE];
- char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
- char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];
- char filerest[FILENAMELENGTH];
- char fileregp[FILENAMELENGTH];
- char path[80],pathc[80],pathcd[80],pathtot[80];
- int firstobs=1, lastobs=10;
- int sdeb, sfin; /* Status at beginning and end */
- int c, h , cpt,l;
- int ju,jl, mi;
- int i1,j1, k1,jk,aa,bb, stepsize;
- int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
-
- int hstepm, nhstepm;
- double bage, fage, age, agelim, agebase;
- double ftolpl=FTOL;
- double **prlim;
- double *severity;
- double ***param; /* Matrix of parameters */
- double *p;
- double **matcov; /* Matrix of covariance */
- double ***delti3; /* Scale */
- double *delti; /* Scale */
- double ***eij, ***vareij;
- double **varpl; /* Variances of prevalence limits by age */
- double *epj, vepp;
- char version[80]="Imach version 0.64, May 2000, INED-EUROREVES ";
- char *alph[]={"a","a","b","c","d","e"}, str[4];
- char z[1]="c";
-#include Computing and drawing one step probabilities with their confidence intervals
Matrix of variance-covariance of pairs of step probabilities (drawings)
Matrix of variance-covariance of pairs of step probabilities
\n\
+ file %s
\n",optionfilehtmcov);
+ fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p
\n");
+ fprintf(fichtmcov,"\n
Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \
+It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \
+would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 \
+standard deviations wide on each axis.
\
+ Now, if both incidences are correlated (usual case) we diagonalised the inverse of the covariance matrix\
+ and made the appropriate rotation to look at the uncorrelated principal directions.
\
+To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.
\n");
+
+ cov[1]=1;
+ tj=cptcoveff;
+ if (cptcovn<1) {tj=1;ncodemax[1]=1;}
+ j1=0;
+ for(t=1; t<=tj;t++){
+ for(i1=1; i1<=ncodemax[t];i1++){
+ j1++;
+ if (cptcovn>0) {
+ fprintf(ficresprob, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficresprob, "**********\n#\n");
+ fprintf(ficresprobcov, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficresprobcov, "**********\n#\n");
+
+ fprintf(ficgp, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficgp, "**********\n#\n");
+
+
+ fprintf(fichtmcov, "\n
********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(fichtmcov, "**********\n
");
+
+ fprintf(ficresprobcor, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficresprobcor, "**********\n#");
+ }
+
+ for (age=bage; age<=fage; age ++){
+ cov[2]=age;
+ for (k=1; k<=cptcovn;k++) {
+ cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
+ }
+ for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
+ for (k=1; k<=cptcovprod;k++)
+ cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+
+ gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
+ trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
+ gp=vector(1,(nlstate)*(nlstate+ndeath));
+ gm=vector(1,(nlstate)*(nlstate+ndeath));
+
+ for(theta=1; theta <=npar; theta++){
+ for(i=1; i<=npar; i++)
+ xp[i] = x[i] + (i==theta ?delti[theta]:(double)0);
+
+ pmij(pmmij,cov,ncovmodel,xp,nlstate);
+
+ k=0;
+ for(i=1; i<= (nlstate); i++){
+ for(j=1; j<=(nlstate+ndeath);j++){
+ k=k+1;
+ gp[k]=pmmij[i][j];
+ }
+ }
+
+ for(i=1; i<=npar; i++)
+ xp[i] = x[i] - (i==theta ?delti[theta]:(double)0);
+
+ pmij(pmmij,cov,ncovmodel,xp,nlstate);
+ k=0;
+ for(i=1; i<=(nlstate); i++){
+ for(j=1; j<=(nlstate+ndeath);j++){
+ k=k+1;
+ gm[k]=pmmij[i][j];
+ }
+ }
+
+ for(i=1; i<= (nlstate)*(nlstate+ndeath); i++)
+ gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta];
+ }
- /* long total_usecs;
- struct timeval start_time, end_time;
-
- gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
+ for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)
+ for(theta=1; theta <=npar; theta++)
+ trgradg[j][theta]=gradg[theta][j];
+
+ matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov);
+ matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
+ free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
+ free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
+ free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+ free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+
+ pmij(pmmij,cov,ncovmodel,x,nlstate);
+
+ k=0;
+ for(i=1; i<=(nlstate); i++){
+ for(j=1; j<=(nlstate+ndeath);j++){
+ k=k+1;
+ mu[k][(int) age]=pmmij[i][j];
+ }
+ }
+ for(i=1;i<=(nlstate)*(nlstate+ndeath);i++)
+ for(j=1;j<=(nlstate)*(nlstate+ndeath);j++)
+ varpij[i][j][(int)age] = doldm[i][j];
+
+ /*printf("\n%d ",(int)age);
+ for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+ printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+ fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+ }*/
+
+ fprintf(ficresprob,"\n%d ",(int)age);
+ fprintf(ficresprobcov,"\n%d ",(int)age);
+ fprintf(ficresprobcor,"\n%d ",(int)age);
+
+ for (i=1; i<=(nlstate)*(nlstate+ndeath);i++)
+ fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age]));
+ for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+ fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]);
+ fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]);
+ }
+ i=0;
+ for (k=1; k<=(nlstate);k++){
+ for (l=1; l<=(nlstate+ndeath);l++){
+ i=i++;
+ fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
+ fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
+ for (j=1; j<=i;j++){
+ fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
+ fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
+ }
+ }
+ }/* end of loop for state */
+ } /* end of loop for age */
+ /* Confidence intervalle of pij */
+ /*
+ fprintf(ficgp,"\nunset parametric;unset label");
+ fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
+ fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+ fprintf(fichtm,"\n
Probability with confidence intervals expressed in year-1 :pijgr%s.png, ",optionfilefiname,optionfilefiname);
+ fprintf(fichtm,"\n
",optionfilefiname);
+ fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
+ fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
+ */
+
+ /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
+ first1=1;
+ for (k2=1; k2<=(nlstate);k2++){
+ for (l2=1; l2<=(nlstate+ndeath);l2++){
+ if(l2==k2) continue;
+ j=(k2-1)*(nlstate+ndeath)+l2;
+ for (k1=1; k1<=(nlstate);k1++){
+ for (l1=1; l1<=(nlstate+ndeath);l1++){
+ if(l1==k1) continue;
+ i=(k1-1)*(nlstate+ndeath)+l1;
+ if(i<=j) continue;
+ for (age=bage; age<=fage; age ++){
+ if ((int)age %5==0){
+ v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;
+ v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM;
+ cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;
+ mu1=mu[i][(int) age]/stepm*YEARM ;
+ mu2=mu[j][(int) age]/stepm*YEARM;
+ c12=cv12/sqrt(v1*v2);
+ /* Computing eigen value of matrix of covariance */
+ lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
+ lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
+ /* Eigen vectors */
+ v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
+ /*v21=sqrt(1.-v11*v11); *//* error */
+ v21=(lc1-v1)/cv12*v11;
+ v12=-v21;
+ v22=v11;
+ tnalp=v21/v11;
+ if(first1==1){
+ first1=0;
+ printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
+ }
+ fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
+ /*printf(fignu*/
+ /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */
+ /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
+ if(first==1){
+ first=0;
+ fprintf(ficgp,"\nset parametric;unset label");
+ fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
+ fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+ fprintf(fichtmcov,"\n
Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\
+ :\
+%s%d%1d%1d-%1d%1d.png, ",k1,l1,k2,l2,\
+ subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\
+ subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ fprintf(fichtmcov,"\n
",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ fprintf(fichtmcov,"\n
Correlation at age %d (%.3f),",(int) age, c12);
+ fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
+ fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
+ fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\
+ mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\
+ mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
+ }else{
+ first=0;
+ fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
+ fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
+ fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
+ fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\
+ mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\
+ mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
+ }/* if first */
+ } /* age mod 5 */
+ } /* end loop age */
+ fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ first=1;
+ } /*l12 */
+ } /* k12 */
+ } /*l1 */
+ }/* k1 */
+ } /* loop covariates */
+ }
+ free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
+ free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
+ free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
+ free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar);
+ free_vector(xp,1,npar);
+ fclose(ficresprob);
+ fclose(ficresprobcov);
+ fclose(ficresprobcor);
+ fflush(ficgp);
+ fflush(fichtmcov);
+}
- printf("\nIMACH, Version 0.64");
- printf("\nEnter the parameter file name: ");
-#define windows 1
-#ifdef windows
- scanf("%s",pathtot);
- getcwd(pathcd, size);
- cut(path,optionfile,pathtot);
- chdir(path);
- replace(pathc,path);
-#endif
-#ifdef unix
- scanf("%s",optionfile);
-#endif
-/*-------- arguments in the command line --------*/
+/******************* 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 rfileres[],\
+ int popforecast, int estepm ,\
+ double jprev1, double mprev1,double anprev1, \
+ double jprev2, double mprev2,double anprev2){
+ int jj1, k1, i1, cpt;
+
+ fprintf(fichtm,"");
+ fprintf(fichtm,"Result files (first order: no variance)
\n \
+ - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s
\n ",
+ jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));
+ fprintf(fichtm,"\
+ - Estimated transition probabilities over %d (stepm) months: %s
\n ",
+ stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
+ fprintf(fichtm,"\
+ - Period (stable) prevalence in each health state: %s
\n",
+ subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
+ fprintf(fichtm,"\
+ - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
+ %s
\n",
+ estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
+ fprintf(fichtm,"\
+ - Population projections by age and states: \
+ %s
\n
");
+
+
+ fprintf(fichtm,"\
+\n
************ Results for covariates");
+ for (cpt=1; cpt<=cptcoveff;cpt++)
+ fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
+ fprintf(fichtm," ************\n
");
+ }
+ /* Pij */
+ fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png
\
+",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ /* Quasi-incidences */
+ fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d2.png
\
+",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ /* Period (stable) prevalence in each health state */
+ for(cpt=1; cpt
\
+",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+ }
+ for(cpt=1; cpt<=nlstate;cpt++) {
+ fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : %s%d%d.png
\
+",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+ }
+ } /* end i1 */
+ }/* End k1 */
+ fprintf(fichtm," Result files (second order: variances)
\n\
+ - Parameter file with estimated parameters and covariance matrix: %s
\n", rfileres,rfileres);
+
+ fprintf(fichtm," - Variance of one-step probabilities: %s
\n",
+ subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
+ fprintf(fichtm,"\
+ - Variance-covariance of one-step probabilities: %s
\n",
+ subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"));
+
+ fprintf(fichtm,"\
+ - Correlation matrix of one-step probabilities: %s
\n",
+ subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
+ fprintf(fichtm,"\
+ - Variances and covariances of health expectancies by age and initial health status (cov(eij,ekl)(estepm=%2d months): \
+ %s
\n
\n",
+ estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
+ fprintf(fichtm,"\
+ - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
\n",
+ estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
+ fprintf(fichtm,"\
+ - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
\n",
+ estepm, subdirf2(fileres,"t"),subdirf2(fileres,"t"));
+ fprintf(fichtm,"\
+ - Standard deviation of period (stable) prevalences: %s
\n",\
+ subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
+
+/* if(popforecast==1) fprintf(fichtm,"\n */
+/* - Prevalences forecasting: f%s
\n */
+/* - Population forecasting (if popforecast=1): pop%s
\n */
+/*
",fileres,fileres,fileres,fileres); */
+/* else */
+/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)
\n",popforecast, stepm, model); */
+ fflush(fichtm);
+ fprintf(fichtm,"
");
+ fflush(fichtm);
+}
- strcpy(fileres,"r");
- strcat(fileres, optionfile);
+/******************* Gnuplot file **************/
+void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
- /*---------arguments file --------*/
+ char dirfileres[132],optfileres[132];
+ int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
+ int ng=0;
+/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
+/* printf("Problem with file %s",optionfilegnuplot); */
+/* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
+/* } */
- if((ficpar=fopen(optionfile,"r"))==NULL) {
- printf("Problem with optionfile %s\n",optionfile);
- goto end;
- }
+ /*#ifdef windows */
+ fprintf(ficgp,"cd \"%s\" \n",pathc);
+ /*#endif */
+ m=pow(2,cptcoveff);
- strcpy(filereso,"o");
- strcat(filereso,fileres);
- if((ficparo=fopen(filereso,"w"))==NULL) {
- printf("Problem with Output resultfile: %s\n", filereso);goto end;
+ strcpy(dirfileres,optionfilefiname);
+ strcpy(optfileres,"vpl");
+ /* 1eme*/
+ for (cpt=1; cpt<= nlstate ; cpt ++) {
+ for (k1=1; k1<= m ; k1 ++) {
+ fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
+ fprintf(ficgp,"\n#set out \"v%s%d%d.png\" \n",optionfilefiname,cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \n\
+set ylabel \"Probability\" \n\
+set ter png small\n\
+set size 0.65,0.65\n\
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
+
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"\" w l 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
+ }
}
-
-/*--------- index.htm --------*/
-
- if((fichtm=fopen("index.htm","w"))==NULL) {
- printf("Problem with index.htm \n");goto end;
+ /*2 eme*/
+
+ for (k1=1; k1<= m ; k1 ++) {
+ fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),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 ++) {
+ k=2*i;
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
+ }
+ if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
+ else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"\" w l 0,");
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
+ }
+ if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
+ else fprintf(ficgp,"\" t\"\" w l 0,");
+ }
}
+
+ /*3eme*/
+
+ for (k1=1; k1<= m ; k1 ++) {
+ for (cpt=1; cpt<= nlstate ; cpt ++) {
+ /* k=2+nlstate*(2*cpt-2); */
+ k=2+(nlstate+1)*(cpt-1);
+ fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
+ fprintf(ficgp,"set ter png small\n\
+set size 0.65,0.65\n\
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),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 ++) {
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+i,cpt,i+1);
+ /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
+
+ }
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+nlstate,cpt);
+ }
+ }
+
+ /* CV preval stable (period) */
+ for (k1=1; k1<= m ; k1 ++) {
+ for (cpt=1; cpt<=nlstate ; cpt ++) {
+ k=3;
+ fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
+set ter png small\nset size 0.65,0.65\n\
+unset log y\n\
+plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,subdirf2(fileres,"pij"),k1,k+cpt+1,k+1);
+
+ for (i=1; i< nlstate ; i ++)
+ fprintf(ficgp,"+$%d",k+i+1);
+ fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
+
+ l=3+(nlstate+ndeath)*cpt;
+ fprintf(ficgp,",\"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",subdirf2(fileres,"pij"),k1,l+cpt+1,l+1);
+ for (i=1; i< nlstate ; i ++) {
+ l=3+(nlstate+ndeath)*cpt;
+ fprintf(ficgp,"+$%d",l+i+1);
+ }
+ fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);
+ }
+ }
+
+ /* proba elementaires */
+ for(i=1,jk=1; i <=nlstate; i++){
+ for(k=1; k <=(nlstate+ndeath); k++){
+ if (k != i) {
+ for(j=1; j <=ncovmodel; j++){
+ fprintf(ficgp,"p%d=%f ",jk,p[jk]);
+ jk++;
+ fprintf(ficgp,"\n");
+ }
+ }
+ }
+ }
- fprintf(fichtm,"
************ Results for covariates");
+ for (cpt=1; cpt<=cptcoveff;cpt++)
+ fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
+ fprintf(fichtm," ************\n
");
+ }
+ for(cpt=1; cpt<=nlstate;cpt++) {
+ fprintf(fichtm,"
- Observed (cross-sectional) and period (incidence based) \
+prevalence (with 95%% confidence interval) in state (%d): %s%d%d.png
\
+",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
+ }
+ fprintf(fichtm,"\n
- Total life expectancy by age and \
+health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
+true period expectancies (those weighted with period prevalences are also\
+ drawn in addition to the population based expectancies computed using\
+ observed and cahotic prevalences: %s%d.png
\
+",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
+ } /* end i1 */
+ }/* End k1 */
+ fprintf(fichtm,"
",strtok(optionfile, "."),strtok(optionfile, "."));
-
+ for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
+ for(jk=1; jk <=m; jk++) {
+ fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng);
+ if (ng==2)
+ fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
+ else
+ fprintf(ficgp,"\nset title \"Probability\"\n");
+ fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar);
+ i=1;
+ for(k2=1; k2<=nlstate; k2++) {
+ k3=i;
+ for(k=1; k<=(nlstate+ndeath); k++) {
+ if (k != k2){
+ if(ng==2)
+ fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+ else
+ fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
+ ij=1;
+ for(j=3; j <=ncovmodel; j++) {
+ if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+ fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ ij++;
+ }
+ else
+ fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+ }
+ fprintf(ficgp,")/(1");
+
+ for(k1=1; k1 <=nlstate; k1++){
+ fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ ij=1;
+ for(j=3; j <=ncovmodel; j++){
+ if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ 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;
+ }
+ } /* end k */
+ } /* end k2 */
+ } /* end jk */
+ } /* end ng */
+ fflush(ficgp);
+} /* end gnuplot */
+
+
+/*************** Moving average **************/
+int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){
+
+ int i, cpt, cptcod;
+ int modcovmax =1;
+ int mobilavrange, mob;
+ double age;
+
+ modcovmax=2*cptcoveff;/* Max number of modalities. We suppose
+ a covariate has 2 modalities */
+ if (cptcovn<1) modcovmax=1; /* At least 1 pass */
+
+ if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
+ if(mobilav==1) mobilavrange=5; /* default */
+ else mobilavrange=mobilav;
+ for (age=bage; age<=fage; age++)
+ for (i=1; i<=nlstate;i++)
+ for (cptcod=1;cptcod<=modcovmax;cptcod++)
+ mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+ /* We keep the original values on the extreme ages bage, fage and for
+ fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
+ we use a 5 terms etc. until the borders are no more concerned.
+ */
+ for (mob=3;mob <=mobilavrange;mob=mob+2){
+ for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
+ for (i=1; i<=nlstate;i++){
+ for (cptcod=1;cptcod<=modcovmax;cptcod++){
+ mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
+ for (cpt=1;cpt<=(mob-1)/2;cpt++){
+ mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
+ mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
+ }
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
+ }
+ }
+ }/* end age */
+ }/* end mob */
+ }else return -1;
+ return 0;
+}/* End movingaverage */
-fclose(fichtm);
- /* Reads comments: lines beginning with '#' */
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
- puts(line);
- fputs(line,ficparo);
+/************** Forecasting ******************/
+prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
+ /* proj1, year, month, day of starting projection
+ agemin, agemax range of age
+ dateprev1 dateprev2 range of dates during which prevalence is computed
+ anproj2 year of en of projection (same day and month as proj1).
+ */
+ int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1;
+ int *popage;
+ double agec; /* generic age */
+ double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double *popeffectif,*popcount;
+ double ***p3mat;
+ double ***mobaverage;
+ char fileresf[FILENAMELENGTH];
+
+ agelim=AGESUP;
+ prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+
+ strcpy(fileresf,"f");
+ strcat(fileresf,fileres);
+ if((ficresf=fopen(fileresf,"w"))==NULL) {
+ printf("Problem with forecast resultfile: %s\n", fileresf);
+ fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
+ }
+ printf("Computing forecasting: result on file '%s' \n", fileresf);
+ fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf);
+
+ if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+ if (mobilav!=0) {
+ mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
+ fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+ printf(" Error in movingaverage mobilav=%d\n",mobilav);
+ }
}
- ungetc(c,ficpar);
- 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\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt);
- 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\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt);
- 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\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt);
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ if (stepm<=12) stepsize=1;
+ if(estepm < stepm){
+ printf ("Problem %d lower than %d\n",estepm, stepm);
+ }
+ else hstepm=estepm;
+
+ hstepm=hstepm/stepm;
+ yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
+ fractional in yp1 */
+ anprojmean=yp;
+ yp2=modf((yp1*12),&yp);
+ mprojmean=yp;
+ yp1=modf((yp2*30.5),&yp);
+ jprojmean=yp;
+ if(jprojmean==0) jprojmean=1;
+ if(mprojmean==0) jprojmean=1;
+
+ i1=cptcoveff;
+ if (cptcovn < 1){i1=1;}
+
+ fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+
+ fprintf(ficresf,"#****** Routine prevforecast **\n");
+
+/* if (h==(int)(YEARM*yearp)){ */
+ for(cptcov=1, k=0;cptcov<=i1;cptcov++){
+ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+ k=k+1;
+ fprintf(ficresf,"\n#******");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ }
+ fprintf(ficresf,"******\n");
+ fprintf(ficresf,"# Covariate valuofcovar yearproj age");
+ for(j=1; j<=nlstate+ndeath;j++){
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresf," p%d%d",i,j);
+ fprintf(ficresf," p.%d",j);
+ }
+ for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
+ fprintf(ficresf,"\n");
+ fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);
+
+ for (agec=fage; agec>=(ageminpar-1); agec--){
+ nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
+ nhstepm = nhstepm/hstepm;
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+
+ for (h=0; h<=nhstepm; h++){
+ if (h*hstepm/YEARM*stepm ==yearp) {
+ fprintf(ficresf,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
+ }
+ for(j=1; j<=nlstate+ndeath;j++) {
+ ppij=0.;
+ for(i=1; i<=nlstate;i++) {
+ if (mobilav==1)
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
+ else {
+ ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
+ }
+ if (h*hstepm/YEARM*stepm== yearp) {
+ fprintf(ficresf," %.3f", p3mat[i][j][h]);
+ }
+ } /* end i */
+ if (h*hstepm/YEARM*stepm==yearp) {
+ fprintf(ficresf," %.3f", ppij);
+ }
+ }/* end j */
+ } /* end h */
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ } /* end agec */
+ } /* end yearp */
+ } /* end cptcod */
+ } /* end cptcov */
+
+ if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+
+ fclose(ficresf);
+}
+
+/************** Forecasting *****not tested NB*************/
+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){
- nvar=ncov-1; /* Suppressing age as a basic covariate */
+ int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
+ int *popage;
+ double calagedatem, agelim, kk1, kk2;
+ double *popeffectif,*popcount;
+ double ***p3mat,***tabpop,***tabpopprev;
+ double ***mobaverage;
+ char filerespop[FILENAMELENGTH];
+
+ tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ agelim=AGESUP;
+ calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
- /* Read guess parameters */
- /* Reads comments: lines beginning with '#' */
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
- puts(line);
- fputs(line,ficparo);
+ prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+
+
+ strcpy(filerespop,"pop");
+ strcat(filerespop,fileres);
+ if((ficrespop=fopen(filerespop,"w"))==NULL) {
+ printf("Problem with forecast resultfile: %s\n", filerespop);
+ fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);
}
- ungetc(c,ficpar);
+ printf("Computing forecasting: result on file '%s' \n", filerespop);
+ fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop);
+
+ if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+ if (mobilav!=0) {
+ mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
+ fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+ printf(" Error in movingaverage mobilav=%d\n",mobilav);
+ }
+ }
+
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ if (stepm<=12) stepsize=1;
- param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncov);
- for(i=1; i <=nlstate; i++)
- for(j=1; j <=nlstate+ndeath-1; j++){
- fscanf(ficpar,"%1d%1d",&i1,&j1);
- fprintf(ficparo,"%1d%1d",i1,j1);
+ agelim=AGESUP;
+
+ hstepm=1;
+ hstepm=hstepm/stepm;
+
+ if (popforecast==1) {
+ if((ficpop=fopen(popfile,"r"))==NULL) {
+ printf("Problem with population file : %s\n",popfile);exit(0);
+ fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0);
+ }
+ popage=ivector(0,AGESUP);
+ popeffectif=vector(0,AGESUP);
+ popcount=vector(0,AGESUP);
+
+ i=1;
+ while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1;
+
+ imx=i;
+ for (i=1; i
\n
- - Observed prevalence in each state: p%s
\n
-- Estimated parameters and the covariance matrix: %s
- - Stationary prevalence in each state: pl%s
- - Transition probabilities: pij%s
- - Copy of the parameter file: o%s
- - 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
",fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
-
- fprintf(fichtm,"
");
-
-for(cpt=1; cpt
",strtok(optionfile, "."),strtok(optionfile, "."),cpt);
- for(cpt=1; cpt<=nlstate;cpt++)
- fprintf(fichtm,"- Observed and stationary prevalence (with confident
-interval) in state (%d): v%s%d.gif
-
",cpt,strtok(optionfile, "."),cpt,strtok(optionfile, "."),cpt);
-
- for(cpt=1; cpt<=nlstate;cpt++)
- fprintf(fichtm,"- Health life expectancies by age and initial health state (%d): exp%s%d.gif
-
",cpt,strtok(optionfile, "."),cpt,strtok(optionfile, "."),cpt);
-
- fprintf(fichtm,"- Total life expectancy by age and
- health expectancies in states (1) and (2): e%s.gif
-
");
+
+fprintf(fichtm,"Result files
\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):
");
+ fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year
",p[1],p[2],agegomp);
+ for (i=1;i<=2;i++)
+ fprintf(fichtm," p[%d] = %lf [%f ; %f]
\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
+ fprintf(fichtm,"
");
+ fprintf(fichtm,"Life table
\n
");
+
+ fprintf(fichtm,"\nAge l
");
+
+ for (k=agegomp;k<(agemortsup-2);k++)
+ fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf
\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);
+
+
+ fflush(fichtm);
+}
+
+/******************* Gnuplot file **************/
+void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+
+ char dirfileres[132],optfileres[132];
+ int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
+ int ng;
+
+
+ /*#ifdef windows */
+ fprintf(ficgp,"cd \"%s\" \n",pathc);
+ /*#endif */
+
+
+ strcpy(dirfileres,optionfilefiname);
+ strcpy(optfileres,"vpl");
+ fprintf(ficgp,"set out \"graphmort.png\"\n ");
+ fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n ");
+ fprintf(ficgp, "set ter png small\n set log y\n");
+ fprintf(ficgp, "set size 0.65,0.65\n");
+ fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp);
+
+}
+
+
+
+
+
+/***********************************************/
+/**************** Main Program *****************/
+/***********************************************/
+
+int main(int argc, char *argv[])
+{
+ int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
+ int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
+ int linei, month, year,iout;
+ int jj, ll, li, lj, lk, imk;
+ int numlinepar=0; /* Current linenumber of parameter file */
+ int itimes;
+ int NDIM=2;
+ int vpopbased=0;
+
+ char ca[32], cb[32], cc[32];
+ char dummy[]=" ";
+ /* FILE *fichtm; *//* Html File */
+ /* FILE *ficgp;*/ /*Gnuplot File */
+ struct stat info;
+ double agedeb, agefin,hf;
+ double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
+
+ double fret;
+ double **xi,tmp,delta;
+
+ double dum; /* Dummy variable */
+ double ***p3mat;
+ double ***mobaverage;
+ int *indx;
+ char line[MAXLINE], linepar[MAXLINE];
+ char linetmp[MAXLINE];
+ char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
+ char pathr[MAXLINE], pathimach[MAXLINE];
+ char **bp, *tok, *val; /* pathtot */
+ int firstobs=1, lastobs=10;
+ int sdeb, sfin; /* Status at beginning and end */
+ 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,*tab;
+ int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
+ int mobilav=0,popforecast=0;
+ int hstepm, nhstepm;
+ int agemortsup;
+ float sumlpop=0.;
+ double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
+ double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;
+
+ double bage, fage, age, agelim, agebase;
+ double ftolpl=FTOL;
+ double **prlim;
+ double *severity;
+ double ***param; /* Matrix of parameters */
+ double *p;
+ double **matcov; /* Matrix of covariance */
+ double ***delti3; /* Scale */
+ double *delti; /* Scale */
+ double ***eij, ***vareij;
+ double **varpl; /* Variances of prevalence limits by age */
+ double *epj, vepp;
+ double kk1, kk2;
+ double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
+ double **ximort;
+ char *alph[]={"a","a","b","c","d","e"}, str[4];
+ int *dcwave;
+
+ char z[1]="c", occ;
+
+ char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
+ char *strt, strtend[80];
+ char *stratrunc;
+ int lstra;
+
+ long total_usecs;
+
+/* setlocale (LC_ALL, ""); */
+/* bindtextdomain (PACKAGE, LOCALEDIR); */
+/* textdomain (PACKAGE); */
+/* setlocale (LC_CTYPE, ""); */
+/* setlocale (LC_MESSAGES, ""); */
+
+ /* gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
+ (void) gettimeofday(&start_time,&tzp);
+ curr_time=start_time;
+ tm = *localtime(&start_time.tv_sec);
+ tmg = *gmtime(&start_time.tv_sec);
+ strcpy(strstart,asctime(&tm));
+
+/* printf("Localtime (at start)=%s",strstart); */
+/* tp.tv_sec = tp.tv_sec +86400; */
+/* tm = *localtime(&start_time.tv_sec); */
+/* tmg.tm_year=tmg.tm_year +dsign*dyear; */
+/* tmg.tm_mon=tmg.tm_mon +dsign*dmonth; */
+/* tmg.tm_hour=tmg.tm_hour + 1; */
+/* tp.tv_sec = mktime(&tmg); */
+/* strt=asctime(&tmg); */
+/* printf("Time(after) =%s",strstart); */
+/* (void) time (&time_value);
+* printf("time=%d,t-=%d\n",time_value,time_value-86400);
+* tm = *localtime(&time_value);
+* strstart=asctime(&tm);
+* printf("tim_value=%d,asctime=%s\n",time_value,strstart);
+*/
+
+ nberr=0; /* Number of errors and warnings */
+ nbwarn=0;
+ getcwd(pathcd, size);
+
+ printf("\n%s\n%s",version,fullversion);
+ if(argc <=1){
+ printf("\nEnter the parameter file name: ");
+ fgets(pathr,FILENAMELENGTH,stdin);
+ i=strlen(pathr);
+ if(pathr[i-1]=='\n')
+ pathr[i-1]='\0';
+ for (tok = pathr; tok != NULL; ){
+ printf("Pathr |%s|\n",pathr);
+ while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
+ printf("val= |%s| pathr=%s\n",val,pathr);
+ strcpy (pathtot, val);
+ if(pathr[0] == '\0') break; /* Dirty */
+ }
+ }
+ else{
+ strcpy(pathtot,argv[1]);
+ }
+ /*if(getcwd(pathcd, MAXLINE)!= NULL)printf ("Error pathcd\n");*/
+ /*cygwin_split_path(pathtot,path,optionfile);
+ printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
+ /* cutv(path,optionfile,pathtot,'\\');*/
+
+ /* Split argv[0], imach program to get pathimach */
+ printf("\nargv[0]=%s argv[1]=%s, \n",argv[0],argv[1]);
+ split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
+ printf("\nargv[0]=%s pathimach=%s, \noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
+ /* strcpy(pathimach,argv[0]); */
+ /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */
+ split(pathtot,path,optionfile,optionfilext,optionfilefiname);
+ printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
+ chdir(path); /* Can be a relative path */
+ if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */
+ printf("Current directory %s!\n",pathcd);
+ strcpy(command,"mkdir ");
+ strcat(command,optionfilefiname);
+ if((outcmd=system(command)) != 0){
+ printf("Problem creating directory or it already exists %s%s, err=%d\n",path,optionfilefiname,outcmd);
+ /* fprintf(ficlog,"Problem creating directory %s%s\n",path,optionfilefiname); */
+ /* fclose(ficlog); */
+/* exit(1); */
+ }
+/* if((imk=mkdir(optionfilefiname))<0){ */
+/* perror("mkdir"); */
+/* } */
+
+ /*-------- arguments in the command line --------*/
+
+ /* Log file */
+ strcat(filelog, optionfilefiname);
+ strcat(filelog,".log"); /* */
+ if((ficlog=fopen(filelog,"w"))==NULL) {
+ printf("Problem with logfile %s\n",filelog);
+ goto end;
+ }
+ fprintf(ficlog,"Log filename:%s\n",filelog);
+ fprintf(ficlog,"\n%s\n%s",version,fullversion);
+ fprintf(ficlog,"\nEnter the parameter file name: \n");
+ fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
+ path=%s \n\
+ optionfile=%s\n\
+ optionfilext=%s\n\
+ optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
+
+ printf("Local time (at start):%s",strstart);
+ fprintf(ficlog,"Local time (at start): %s",strstart);
+ fflush(ficlog);
+/* (void) gettimeofday(&curr_time,&tzp); */
+/* printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */
+
+ /* */
+ strcpy(fileres,"r");
+ strcat(fileres, optionfilefiname);
+ strcat(fileres,".txt"); /* Other files have txt extension */
+
+ /*---------arguments file --------*/
+
+ if((ficpar=fopen(optionfile,"r"))==NULL) {
+ printf("Problem with optionfile %s\n",optionfile);
+ fprintf(ficlog,"Problem with optionfile %s\n",optionfile);
+ fflush(ficlog);
+ goto end;
+ }
+
+
+
+ strcpy(filereso,"o");
+ strcat(filereso,fileres);
+ if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */
+ printf("Problem with Output resultfile: %s\n", filereso);
+ fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);
+ fflush(ficlog);
+ goto end;
+ }
+
/* Reads comments: lines beginning with '#' */
+ numlinepar=0;
while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);
+ numlinepar++;
puts(line);
fputs(line,ficparo);
+ fputs(line,ficlog);
}
ungetc(c,ficpar);
-
- matcov=matrix(1,npar,1,npar);
- for(i=1; i <=npar; i++){
- fscanf(ficpar,"%s",&str);
- printf("%s",str);
- fprintf(ficparo,"%s",str);
- for(j=1; j <=i; j++){
- fscanf(ficpar," %le",&matcov[i][j]);
- printf(" %.5le",matcov[i][j]);
- fprintf(ficparo," %.5le",matcov[i][j]);
- }
- fscanf(ficpar,"\n");
- printf("\n");
- fprintf(ficparo,"\n");
+
+ fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
+ numlinepar++;
+ printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+ fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+ fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+ fflush(ficlog);
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ numlinepar++;
+ puts(line);
+ fputs(line,ficparo);
+ fputs(line,ficlog);
}
- for(i=1; i <=npar; i++)
- for(j=i+1;j<=npar;j++)
- matcov[i][j]=matcov[j][i];
+ ungetc(c,ficpar);
+
- printf("\n");
+ covar=matrix(0,NCOVMAX,1,n);
+ cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/
+ if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
+ /* where is ncovprod ?*/
+ ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/
+ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
+ nforces= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */
+ npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters like aij*/
+ if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){
+ printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX);
+ fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX);
+ fflush(stdout);
+ fclose (ficlog);
+ goto end;
+ }
+ delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
+ delti=delti3[1][1];
+ /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
+ if(mle==-1){ /* Print a wizard for help writing covariance matrix */
+ prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
+ printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
+ fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
+ free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
+ fclose (ficparo);
+ fclose (ficlog);
+ goto end;
+ exit(0);
+ }
+ else if(mle==-3) {
+ prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
+ printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
+ fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
+ param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
+ matcov=matrix(1,npar,1,npar);
+ }
+ else{
+ /* Read guess parameters */
+ /* Reads comments: lines beginning with '#' */
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ numlinepar++;
+ puts(line);
+ fputs(line,ficparo);
+ fputs(line,ficlog);
+ }
+ ungetc(c,ficpar);
+
+ param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
+ for(i=1; i <=nlstate; i++){
+ j=0;
+ for(jj=1; jj <=nlstate+ndeath; jj++){
+ if(jj==i) continue;
+ j++;
+ fscanf(ficpar,"%1d%1d",&i1,&j1);
+ if ((i1 != i) && (j1 != j)){
+ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
+It might be a problem of design; if ncovcol and the model are correct\n \
+run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
+ exit(1);
+ }
+ fprintf(ficparo,"%1d%1d",i1,j1);
+ if(mle==1)
+ printf("%1d%1d",i,j);
+ fprintf(ficlog,"%1d%1d",i,j);
+ for(k=1; k<=ncovmodel;k++){
+ fscanf(ficpar," %lf",¶m[i][j][k]);
+ if(mle==1){
+ printf(" %lf",param[i][j][k]);
+ fprintf(ficlog," %lf",param[i][j][k]);
+ }
+ else
+ fprintf(ficlog," %lf",param[i][j][k]);
+ fprintf(ficparo," %lf",param[i][j][k]);
+ }
+ fscanf(ficpar,"\n");
+ numlinepar++;
+ if(mle==1)
+ printf("\n");
+ fprintf(ficlog,"\n");
+ fprintf(ficparo,"\n");
+ }
+ }
+ fflush(ficlog);
+
+ p=param[1][1];
+
+ /* Reads comments: lines beginning with '#' */
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ numlinepar++;
+ puts(line);
+ fputs(line,ficparo);
+ fputs(line,ficlog);
+ }
+ ungetc(c,ficpar);
+
+ for(i=1; i <=nlstate; i++){
+ for(j=1; j <=nlstate+ndeath-1; j++){
+ fscanf(ficpar,"%1d%1d",&i1,&j1);
+ if ((i1-i)*(j1-j)!=0){
+ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
+ exit(1);
+ }
+ printf("%1d%1d",i,j);
+ fprintf(ficparo,"%1d%1d",i1,j1);
+ fprintf(ficlog,"%1d%1d",i1,j1);
+ for(k=1; k<=ncovmodel;k++){
+ fscanf(ficpar,"%le",&delti3[i][j][k]);
+ printf(" %le",delti3[i][j][k]);
+ fprintf(ficparo," %le",delti3[i][j][k]);
+ fprintf(ficlog," %le",delti3[i][j][k]);
+ }
+ fscanf(ficpar,"\n");
+ numlinepar++;
+ printf("\n");
+ fprintf(ficparo,"\n");
+ fprintf(ficlog,"\n");
+ }
+ }
+ fflush(ficlog);
+
+ delti=delti3[1][1];
+
+
+ /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */
+ /* Reads comments: lines beginning with '#' */
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ numlinepar++;
+ puts(line);
+ fputs(line,ficparo);
+ fputs(line,ficlog);
+ }
+ ungetc(c,ficpar);
- if(mle==1){
- /*-------- data file ----------*/
- if((ficres =fopen(fileres,"w"))==NULL) {
- printf("Problem with resultfile: %s\n", fileres);goto end;
+ matcov=matrix(1,npar,1,npar);
+ for(i=1; i <=npar; i++)
+ for(j=1; j <=npar; j++) matcov[i][j]=0.;
+
+ for(i=1; i <=npar; i++){
+ fscanf(ficpar,"%s",&str);
+ if(mle==1)
+ printf("%s",str);
+ fprintf(ficlog,"%s",str);
+ fprintf(ficparo,"%s",str);
+ for(j=1; j <=i; j++){
+ fscanf(ficpar," %le",&matcov[i][j]);
+ if(mle==1){
+ printf(" %.5le",matcov[i][j]);
+ }
+ fprintf(ficlog," %.5le",matcov[i][j]);
+ fprintf(ficparo," %.5le",matcov[i][j]);
+ }
+ fscanf(ficpar,"\n");
+ numlinepar++;
+ if(mle==1)
+ printf("\n");
+ fprintf(ficlog,"\n");
+ fprintf(ficparo,"\n");
+ }
+ for(i=1; i <=npar; i++)
+ for(j=i+1;j<=npar;j++)
+ matcov[i][j]=matcov[j][i];
+
+ if(mle==1)
+ printf("\n");
+ fprintf(ficlog,"\n");
+
+ fflush(ficlog);
+
+ /*-------- Rewriting parameter file ----------*/
+ strcpy(rfileres,"r"); /* "Rparameterfile */
+ strcat(rfileres,optionfilefiname); /* Parameter file first name*/
+ strcat(rfileres,"."); /* */
+ strcat(rfileres,optionfilext); /* Other files have txt extension */
+ if((ficres =fopen(rfileres,"w"))==NULL) {
+ printf("Problem writing new parameter file: %s\n", fileres);goto end;
+ fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
}
fprintf(ficres,"#%s\n",version);
+ } /* End of mle != -3 */
+
+ /*-------- data file ----------*/
+ if((fic=fopen(datafile,"r"))==NULL) {
+ printf("Problem while opening datafile: %s\n", datafile);goto end;
+ fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end;
+ }
+
+ n= lastobs;
+ severity = vector(1,maxwav);
+ outcome=imatrix(1,maxwav+1,1,n);
+ num=lvector(1,n);
+ moisnais=vector(1,n);
+ annais=vector(1,n);
+ moisdc=vector(1,n);
+ andc=vector(1,n);
+ agedc=vector(1,n);
+ cod=ivector(1,n);
+ weight=vector(1,n);
+ for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
+ mint=matrix(1,maxwav,1,n);
+ anint=matrix(1,maxwav,1,n);
+ s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */
+ tab=ivector(1,NCOVMAX);
+ ncodemax=ivector(1,8);
+
+ i=1;
+ linei=0;
+ while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
+ linei=linei+1;
+ for(j=strlen(line); j>=0;j--){ /* Untabifies line */
+ if(line[j] == '\t')
+ line[j] = ' ';
+ }
+ for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){
+ ;
+ };
+ line[j+1]=0; /* Trims blanks at end of line */
+ if(line[0]=='#'){
+ fprintf(ficlog,"Comment line\n%s\n",line);
+ printf("Comment line\n%s\n",line);
+ continue;
+ }
+ trimbb(linetmp,line); /* Trims multiple blanks in line */
+ for (j=0; line[j]!='\0';j++){
+ line[j]=linetmp[j];
+ }
+
+
+ for (j=maxwav;j>=1;j--){
+ cutv(stra, strb,line,' ');
+ if(strb[0]=='.') { /* Missing status */
+ lval=-1;
+ }else{
+ errno=0;
+ lval=strtol(strb,&endptr,10);
+ /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+ goto end;
+ }
+ }
+ s[j][i]=lval;
+
+ strcpy(line,stra);
+ cutv(stra, strb,line,' ');
+ if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+ }
+ else if(iout=sscanf(strb,"%s.") != 0){
+ month=99;
+ year=9999;
+ }else{
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
+ goto end;
+ }
+ anint[j][i]= (double) year;
+ mint[j][i]= (double)month;
+ strcpy(line,stra);
+ } /* ENd Waves */
+
+ cutv(stra, strb,line,' ');
+ if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+ }
+ else if(iout=sscanf(strb,"%s.",dummy) != 0){
+ month=99;
+ year=9999;
+ }else{
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog);
+ goto end;
+ }
+ andc[i]=(double) year;
+ moisdc[i]=(double) month;
+ strcpy(line,stra);
- if((fic=fopen(datafile,"r"))==NULL) {
- printf("Problem with datafile: %s\n", datafile);goto end;
+ cutv(stra, strb,line,' ');
+ if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+ }
+ else if(iout=sscanf(strb,"%s.") != 0){
+ month=99;
+ year=9999;
+ }else{
+ printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);fflush(ficlog);
+ goto end;
}
+ annais[i]=(double)(year);
+ moisnais[i]=(double)(month);
+ strcpy(line,stra);
- n= lastobs;
- severity = vector(1,maxwav);
- outcome=imatrix(1,maxwav+1,1,n);
- num=ivector(1,n);
- moisnais=vector(1,n);
- annais=vector(1,n);
- moisdc=vector(1,n);
- andc=vector(1,n);
- agedc=vector(1,n);
- cod=ivector(1,n);
- weight=vector(1,n);
- for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
- mint=matrix(1,maxwav,1,n);
- anint=matrix(1,maxwav,1,n);
- covar=matrix(1,NCOVMAX,1,n);
- s=imatrix(1,maxwav+1,1,n);
- adl=imatrix(1,maxwav+1,1,n);
- tab=ivector(1,NCOVMAX);
- i=1;
- while (fgets(line, MAXLINE, fic) != NULL) {
- if ((i >= firstobs) && (i
%s \
+
\n\
+Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n",\
+ optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+ }
+
+ fprintf(fichtm,"\n
%s \
+
\n\
+Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n\
+\n\
+
\
+
\n",\
+ optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
+ optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
+ fileres,fileres,\
+ filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
+ fflush(fichtm);
+
+ strcpy(pathr,path);
+ strcat(pathr,optionfilefiname);
+ chdir(optionfilefiname); /* Move to directory named optionfile */
+
+ /* Calculates basic frequencies. Computes observed prevalence at single age
+ and prints on file fileres'p'. */
+ freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart);
+
+ fprintf(fichtm,"\n");
+ fprintf(fichtm,"Parameter files
\n\
+ - Parameter file: %s.%s
\n\
+ - Copy of the parameter file: o%s
\n\
+ - Log file of the run: %s
\n\
+ - Gnuplot file name: %s
\n\
+ - Date and time at start: %s
Total number of observations=%d
\n\
+Youngest age at first (selected) pass %.2f, oldest age %.2f
\n\
+Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\
+ imx,agemin,agemax,jmin,jmax,jmean);
+ pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
- /* For Powell, parameters are in a vector p[] starting at p[1]
- so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
- p=param[1][1]; /* *(*(*(param +1)+1)+0) */
+
+ /* For Powell, parameters are in a vector p[] starting at p[1]
+ so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
+ p=param[1][1]; /* *(*(*(param +1)+1)+0) */
+
+ globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
+
+ if (mle==-3){
+ ximort=matrix(1,NDIM,1,NDIM);
+ cens=ivector(1,n);
+ ageexmed=vector(1,n);
+ agecens=vector(1,n);
+ dcwave=ivector(1,n);
+
+ for (i=1; i<=imx; i++){
+ dcwave[i]=-1;
+ for (m=firstpass; m<=lastpass; m++)
+ if (s[m][i]>nlstate) {
+ dcwave[i]=m;
+ /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
+ break;
+ }
+ }
+
+ for (i=1; i<=imx; i++) {
+ if (wav[i]>0){
+ ageexmed[i]=agev[mw[1][i]][i];
+ j=wav[i];
+ agecens[i]=1.;
+
+ if (ageexmed[i]> 1 && wav[i] > 0){
+ agecens[i]=agev[mw[j][i]][i];
+ cens[i]= 1;
+ }else if (ageexmed[i]< 1)
+ cens[i]= -1;
+ if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
+ cens[i]=0 ;
+ }
+ else cens[i]=-1;
+ }
+
+ for (i=1;i<=NDIM;i++) {
+ for (j=1;j<=NDIM;j++)
+ ximort[i][j]=(i == j ? 1.0 : 0.0);
+ }
+
+ p[1]=0.0268; p[NDIM]=0.083;
+ /*printf("%lf %lf", p[1], p[2]);*/
- mlikeli(ficres,p, npar, ncov, nlstate, ftol, func);
+
+ printf("Powell\n"); fprintf(ficlog,"Powell\n");
+ strcpy(filerespow,"pow-mort");
+ strcat(filerespow,fileres);
+ if((ficrespow=fopen(filerespow,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", filerespow);
+ fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
+ }
+ fprintf(ficrespow,"# Powell\n# iter -2*LL");
+ /* for (i=1;i<=nlstate;i++)
+ for(j=1;j<=nlstate+ndeath;j++)
+ if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
+ */
+ fprintf(ficrespow,"\n");
+
+ powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
+ fclose(ficrespow);
+
+ hesscov(matcov, p, NDIM, delti, 1e-4, gompertz);
+ for(i=1; i <=NDIM; i++)
+ for(j=i+1;j<=NDIM;j++)
+ matcov[i][j]=matcov[j][i];
+
+ printf("\nCovariance matrix\n ");
+ for(i=1; i <=NDIM; i++) {
+ for(j=1;j<=NDIM;j++){
+ printf("%f ",matcov[i][j]);
+ }
+ printf("\n ");
+ }
+
+ printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
+ for (i=1;i<=NDIM;i++)
+ printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
+
+ lsurv=vector(1,AGESUP);
+ lpop=vector(1,AGESUP);
+ tpop=vector(1,AGESUP);
+ lsurv[agegomp]=100000;
+
+ for (k=agegomp;k<=AGESUP;k++) {
+ agemortsup=k;
+ if (p[1]*exp(p[2]*(k-agegomp))>1) break;
+ }
+
+ for (k=agegomp;k
Local time at start %s
Local time at end %s
\n",strstart, strtend);
+ fclose(fichtm);
+ fprintf(fichtmcov,"
Local time at start %s
Local time at end %s
\n",strstart, strtend);
+ fclose(fichtmcov);
+ fclose(ficgp);
+ fclose(ficlog);
/*------ End -----------*/
- end:
-#ifdef windows
- chdir(pathcd);
-#endif
- system("gnuplot graph.gp");
-#ifdef windows
+ printf("Before Current directory %s!\n",pathcd);
+ if(chdir(pathcd) != 0)
+ printf("Can't move to directory %s!\n",path);
+ if(getcwd(pathcd,MAXLINE) > 0)
+ printf("Current directory %s!\n",pathcd);
+ /*strcat(plotcmd,CHARSEPARATOR);*/
+ sprintf(plotcmd,"gnuplot");
+#ifndef UNIX
+ sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);
+#endif
+ if(!stat(plotcmd,&info)){
+ printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
+ if(!stat(getenv("GNUPLOTBIN"),&info)){
+ printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);
+ }else
+ strcpy(pplotcmd,plotcmd);
+#ifdef UNIX
+ strcpy(plotcmd,GNUPLOTPROGRAM);
+ if(!stat(plotcmd,&info)){
+ printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
+ }else
+ strcpy(pplotcmd,plotcmd);
+#endif
+ }else
+ strcpy(pplotcmd,plotcmd);
+
+ sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot);
+ printf("Starting graphs with: %s\n",plotcmd);fflush(stdout);
+
+ if((outcmd=system(plotcmd)) != 0){
+ printf("\n Problem with gnuplot\n");
+ }
+ printf(" Wait...");
while (z[0] != 'q') {
- chdir(pathcd);
- printf("\nType e to edit output files, c to start again, and q for exiting: ");
+ /* chdir(path); */
+ printf("\nType e to edit output files, g to graph again and q for exiting: ");
scanf("%s",z);
- if (z[0] == 'c') system("./imach");
- else if (z[0] == 'e') {
- chdir(path);
- system("index.htm");
+/* if (z[0] == 'c') system("./imach"); */
+ if (z[0] == 'e') {
+ printf("Starting browser with: %s",optionfilehtm);fflush(stdout);
+ system(optionfilehtm);
}
+ else if (z[0] == 'g') system(plotcmd);
else if (z[0] == 'q') exit(0);
}
-#endif
+ end:
+ while (z[0] != 'q') {
+ printf("\nType q for exiting: ");
+ scanf("%s",z);
+ }
}
+
+