--- imach/src/imach.c 2015/09/22 19:45:16 1.202 +++ imach/src/imach.c 2015/12/23 17:18:31 1.217 @@ -1,6 +1,69 @@ -/* $Id: imach.c,v 1.202 2015/09/22 19:45:16 brouard Exp $ +/* $Id: imach.c,v 1.217 2015/12/23 17:18:31 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.217 2015/12/23 17:18:31 brouard + Summary: Experimental backcast + + Revision 1.216 2015/12/18 17:32:11 brouard + Summary: 0.98r4 Warning and status=-2 + + Version 0.98r4 is now: + - displaying an error when status is -1, date of interview unknown and date of death known; + - permitting a status -2 when the vital status is unknown at a known date of right truncation. + Older changes concerning s=-2, dating from 2005 have been supersed. + + Revision 1.215 2015/12/16 08:52:24 brouard + Summary: 0.98r4 working + + Revision 1.214 2015/12/16 06:57:54 brouard + Summary: temporary not working + + Revision 1.213 2015/12/11 18:22:17 brouard + Summary: 0.98r4 + + Revision 1.212 2015/11/21 12:47:24 brouard + Summary: minor typo + + Revision 1.211 2015/11/21 12:41:11 brouard + Summary: 0.98r3 with some graph of projected cross-sectional + + Author: Nicolas Brouard + + Revision 1.210 2015/11/18 17:41:20 brouard + Summary: Start working on projected prevalences + + Revision 1.209 2015/11/17 22:12:03 brouard + Summary: Adding ftolpl parameter + Author: N Brouard + + We had difficulties to get smoothed confidence intervals. It was due + to the period prevalence which wasn't computed accurately. The inner + parameter ftolpl is now an outer parameter of the .imach parameter + file after estepm. If ftolpl is small 1.e-4 and estepm too, + computation are long. + + Revision 1.208 2015/11/17 14:31:57 brouard + Summary: temporary + + Revision 1.207 2015/10/27 17:36:57 brouard + *** empty log message *** + + Revision 1.206 2015/10/24 07:14:11 brouard + *** empty log message *** + + Revision 1.205 2015/10/23 15:50:53 brouard + Summary: 0.98r3 some clarification for graphs on likelihood contributions + + Revision 1.204 2015/10/01 16:20:26 brouard + Summary: Some new graphs of contribution to likelihood + + Revision 1.203 2015/09/30 17:45:14 brouard + Summary: looking at better estimation of the hessian + + Also a better criteria for convergence to the period prevalence And + therefore adding the number of years needed to converge. (The + prevalence in any alive state shold sum to one + Revision 1.202 2015/09/22 19:45:16 brouard Summary: Adding some overall graph on contribution to likelihood. Might change @@ -647,7 +710,10 @@ /* #define DEBUG */ /* #define DEBUGBRENT */ -#define DEBUGLINMIN +/* #define DEBUGLINMIN */ +/* #define DEBUGHESS */ +#define DEBUGHESSIJ +/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */ #define POWELL /* Instead of NLOPT */ #define POWELLF1F3 /* Skip test */ /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ @@ -719,6 +785,8 @@ typedef struct { #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 +/*#define decodtabm(h,k,cptcoveff)= (h <= (1<> (k-1)) & 1) +1 : -1)*/ +#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 #define MAXN 20000 #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 @@ -735,12 +803,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.202 2015/09/22 19:45:16 brouard Exp $ */ +/* $Id: imach.c,v 1.217 2015/12/23 17:18:31 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.202 $ $Date: 2015/09/22 19:45:16 $"; +char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; +char fullversion[]="$Revision: 1.217 $ $Date: 2015/12/23 17:18:31 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -776,7 +844,7 @@ double **matprod2(); /* test */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ /*FILE *fic ; */ /* Used in readdata only */ -FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; +FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficresplb,*ficrespij, *ficrespijb, *ficrest,*ficresf, *ficresfb,*ficrespop; FILE *ficlog, *ficrespow; int globpr=0; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ @@ -799,13 +867,13 @@ char fileresv[FILENAMELENGTH]; FILE *ficresvpl; char fileresvpl[FILENAMELENGTH]; char title[MAXLINE]; -char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; +char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH]; char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH]; char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; char command[FILENAMELENGTH]; int outcmd=0; -char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; +char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filerespijb[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; char fileresu[FILENAMELENGTH]; /* fileres without r in front */ char filelog[FILENAMELENGTH]; /* Log file */ char filerest[FILENAMELENGTH]; @@ -930,7 +998,7 @@ static int split( char *path, char *dirc } /* got dirc from getcwd*/ printf(" DIRC = %s \n",dirc); - } else { /* strip direcotry from path */ + } else { /* strip directory from path */ ss++; /* after this, the filename */ l2 = strlen( ss ); /* length of filename */ if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); @@ -1338,7 +1406,30 @@ char *subdirf3(char fileres[], char *pre strcat(tmpout,fileres); return tmpout; } + +/*************** function subdirfext ***********/ +char *subdirfext(char fileres[], char *preop, char *postop) +{ + + strcpy(tmpout,preop); + strcat(tmpout,fileres); + strcat(tmpout,postop); + return tmpout; +} +/*************** function subdirfext3 ***********/ +char *subdirfext3(char fileres[], char *preop, char *postop) +{ + + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); + strcat(tmpout,preop); + strcat(tmpout,fileres); + strcat(tmpout,postop); + return tmpout; +} + char *asc_diff_time(long time_sec, char ascdiff[]) { long sec_left, days, hours, minutes; @@ -1602,8 +1693,11 @@ void linmin(double p[], double xi[], int double xx,xmin,bx,ax; double fx,fb,fa; - double scale=10., axs, xxs, xxss; /* Scale added for infinity */ - +#ifdef LINMINORIGINAL +#else + double scale=10., axs, xxs; /* Scale added for infinity */ +#endif + ncom=n; pcom=vector(1,n); xicom=vector(1,n); @@ -1613,12 +1707,15 @@ void linmin(double p[], double xi[], int xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */ } - /* axs=0.0; */ - /* xxss=1; /\* 1 and using scale *\/ */ - xxs=1; - /* do{ */ - ax=0.; +#ifdef LINMINORIGINAL + xx=1.; +#else + axs=0.0; + xxs=1.; + do{ xx= xxs; +#endif + ax=0.; mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */ /* brackets with inputs ax=0 and xx=1, but points, pcom=p, and directions values, xicom=xi, are sent via f1dim(x) */ /* xt[x,j]=pcom[j]+x*xicom[j] f(ax) = f(xt(a,j=1,n)) = f(p(j) + 0 * xi(j)) and f(xx) = f(xt(x, j=1,n)) = f(p(j) + 1 * xi(j)) */ @@ -1626,12 +1723,19 @@ void linmin(double p[], double xi[], int /* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */ /* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */ /* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/ - /* if (fx != fx){ */ - /* xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */ - /* printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); */ - /* } */ - /* }while(fx != fx); */ - +#ifdef LINMINORIGINAL +#else + if (fx != fx){ + xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */ + printf("|"); + fprintf(ficlog,"|"); +#ifdef DEBUGLINMIN + printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); +#endif + } + }while(fx != fx); +#endif + #ifdef DEBUGLINMIN printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); @@ -1650,14 +1754,23 @@ void linmin(double p[], double xi[], int fprintf(ficlog,"linmin end "); #endif for (j=1;j<=n;j++) { - /* printf(" before xi[%d]=%12.8f", j,xi[j]); */ - xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ - /* if(xxs <1.0) */ - /* printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); */ +#ifdef LINMINORIGINAL + xi[j] *= xmin; +#else +#ifdef DEBUGLINMIN + if(xxs <1.0) + printf(" before xi[%d]=%12.8f", j,xi[j]); +#endif + xi[j] *= xmin*xxs; /* xi rescaled by xmin and number of loops: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ +#ifdef DEBUGLINMIN + if(xxs <1.0) + printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); +#endif +#endif p[j] += xi[j]; /* Parameters values are updated accordingly */ } - /* printf("\n"); */ #ifdef DEBUGLINMIN + printf("\n"); printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); fprintf(ficlog,"Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); for (j=1;j<=n;j++) { @@ -1668,6 +1781,7 @@ void linmin(double p[], double xi[], int fprintf(ficlog,"\n"); } } +#else #endif free_vector(xicom,1,n); free_vector(pcom,1,n); @@ -1745,10 +1859,10 @@ void powell(double p[], double **xi, int for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */ fptt=(*fret); #ifdef DEBUG - printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); - fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); + printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); + fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); #endif - printf("%d",i);fflush(stdout); /* print direction (parameter) i */ + printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/ /* Outputs are fret(new point p) p is updated and xit rescaled */ @@ -1917,19 +2031,40 @@ void powell(double p[], double **xi, int /**** Prevalence limit (stable or period prevalence) ****************/ -double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) +double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, 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 */ - + matrix by transitions matrix until convergence is reached with precision ftolpl */ + /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ + /* Wx is row vector: population in state 1, population in state 2, population dead */ + /* or prevalence in state 1, prevalence in state 2, 0 */ + /* newm is the matrix after multiplications, its rows are identical at a factor */ + /* Initial matrix pimij */ + /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ + /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ + /* 0, 0 , 1} */ + /* + * and after some iteration: */ + /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */ + /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */ + /* 0, 0 , 1} */ + /* And prevalence by suppressing the deaths are close to identical rows in prlim: */ + /* {0.51571254859325999, 0.4842874514067399, */ + /* 0.51326036147820708, 0.48673963852179264} */ + /* If we start from prlim again, prlim tends to a constant matrix */ + int i, ii,j,k; - double min, max, maxmin, maxmax,sumnew=0.; + double *min, *max, *meandiff, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ double **out, cov[NCOVMAX+1], **pmij(); double **newm; - double agefin, delaymax=100 ; /* Max number of years to converge */ - long int ncvyear=0, ncvloop=0; + double agefin, delaymax=200. ; /* 100 Max number of years to converge */ + int ncvloop=0; + min=vector(1,nlstate); + max=vector(1,nlstate); + meandiff=vector(1,nlstate); + for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); @@ -1967,32 +2102,165 @@ double **prevalim(double **prlim, int nl savm=oldm; oldm=newm; - maxmax=0.; - for(j=1;j<=nlstate;j++){ - min=1.; - max=0.; - for(i=1; i<=nlstate; i++) { - sumnew=0; - for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; + + for(j=1; j<=nlstate; j++){ + max[j]=0.; + min[j]=1.; + } + for(i=1;i<=nlstate;i++){ + sumnew=0; + for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; + for(j=1; j<=nlstate; j++){ prlim[i][j]= newm[i][j]/(1-sumnew); - max=FMAX(max,prlim[i][j]); - min=FMIN(min,prlim[i][j]); - /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */ + max[j]=FMAX(max[j],prlim[i][j]); + min[j]=FMIN(min[j],prlim[i][j]); } - maxmin=max-min; - maxmax=FMAX(maxmax,maxmin); + } + + maxmax=0.; + for(j=1; j<=nlstate; j++){ + meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */ + maxmax=FMAX(maxmax,meandiff[j]); + /* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */ } /* j loop */ + *ncvyear= (int)age- (int)agefin; + /* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ if(maxmax < ftolpl){ - /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age-(int)agefin); */ + /* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */ + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); return prlim; } } /* age loop */ - printf("Warning: the stable prevalence did not converge with the required precision ftolpl=6*10^5*ftol=%g. \n\ -Earliest age to start was %d-%d=%d, ncvloop=%ld, ncvyear=%d\n\ -Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); + /* After some age loop it doesn't converge */ + printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ +Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); + return prlim; /* should not reach here */ } + + /**** Back Prevalence limit (stable or period prevalence) ****************/ + +double **bprevalim(double **bprlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, 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 with precision ftolpl */ + /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ + /* Wx is row vector: population in state 1, population in state 2, population dead */ + /* or prevalence in state 1, prevalence in state 2, 0 */ + /* newm is the matrix after multiplications, its rows are identical at a factor */ + /* Initial matrix pimij */ + /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ + /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ + /* 0, 0 , 1} */ + /* + * and after some iteration: */ + /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */ + /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */ + /* 0, 0 , 1} */ + /* And prevalence by suppressing the deaths are close to identical rows in prlim: */ + /* {0.51571254859325999, 0.4842874514067399, */ + /* 0.51326036147820708, 0.48673963852179264} */ + /* If we start from prlim again, prlim tends to a constant matrix */ + + int i, ii,j,k; + double *min, *max, *meandiff, maxmax,sumnew=0.; + /* double **matprod2(); */ /* test */ + double **out, cov[NCOVMAX+1], **bmij(); + double **newm; + double agefin, delaymax=200. ; /* 100 Max number of years to converge */ + int ncvloop=0; + + min=vector(1,nlstate); + max=vector(1,nlstate); + meandiff=vector(1,nlstate); + + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + } + + cov[1]=1.; + + /* Even if hstepm = 1, at least one multiplication by the unit matrix */ + /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */ + for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ + ncvloop++; + newm=savm; + /* Covariates have to be included here again */ + cov[2]=agefin; + if(nagesqr==1) + cov[3]= agefin*agefin;; + for (k=1; k<=cptcovn;k++) { + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; + /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */ + } + /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */ + for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; + for (k=1; k<=cptcovprod;k++) /* Useless */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; + + /*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]);*/ + /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ + /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */ + out=matprod2(newm, oldm ,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate)); /* Bug Valgrind */ + + savm=oldm; + oldm=newm; + + for(j=1; j<=nlstate; j++){ + max[j]=0.; + min[j]=1.; + } + /* sumnew=0; */ + /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */ + for(j=1; j<=nlstate; j++){ + for(i=1;i<=nlstate;i++){ + /* bprlim[i][j]= newm[i][j]/(1-sumnew); */ + bprlim[i][j]= newm[i][j]; + max[i]=FMAX(max[i],bprlim[i][j]); + min[i]=FMIN(min[i],bprlim[i][j]); + } + } + + maxmax=0.; + for(i=1; i<=nlstate; i++){ + meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */ + maxmax=FMAX(maxmax,meandiff[i]); + /* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */ + } /* j loop */ + *ncvyear= -( (int)age- (int)agefin); + /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */ + if(maxmax < ftolpl){ + printf("OK Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); + return bprlim; + } + } /* age loop */ + /* After some age loop it doesn't converge */ + printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ +Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); + + return bprlim; /* should not reach here */ +} + /*************** transition probabilities ***************/ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) @@ -2075,6 +2343,98 @@ double **pmij(double **ps, double *cov, return ps; } +/*************** transition probabilities ***************/ + +double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) +{ + /* According to parameters values stored in x and the covariate's values stored in cov, + computes the probability to be observed in state j being in state i by appying the + model to the ncovmodel covariates (including constant and age). + lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc] + and, according on how parameters are entered, the position of the coefficient xij(nc) of the + ncth covariate in the global vector x is given by the formula: + j=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel + Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation, + sums on j different of i to get 1-pii/pii, deduces pii, and then all pij. + Outputs ps[i][j] the probability to be observed in j being in j according to + the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij] + */ + double s1, lnpijopii; + /*double t34;*/ + int i,j, nc, ii, jj; + + for(i=1; i<= nlstate; i++){ + for(j=1; ji s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ + } + ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ + } + } + + for(i=1; i<= nlstate; i++){ + s1=0; + for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ + ps[i][i]=1./(s1+1.); + /* Computing other pijs */ + for(j=1; jfunction)(xt); /* p xt[1]@8 is fine */ - /* fret=(*func)(xt); /\* p xt[1]@8 is fine *\/ */ - printf("Function = %.12lf ",fret); - for (j=1;j<=n;j++) printf(" %d %.8lf", j, xt[j]); - printf("\n"); - free_vector(xt,1,n); - return fret; -} -#endif +/************* Higher Back Matrix Product ***************/ -/*************** log-likelihood *************/ -double func( double *x) +double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij ) { - int i, ii, j, k, mi, d, kk; - double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; - double **out; - double sw; /* Sum of weights */ - double lli; /* Individual log likelihood */ - int s1, s2; - double bbh, survp; - long ipmx; + /* Computes the transition matrix starting at age 'age' over + 'nhstepm*hstepm*stepm' months (i.e. until + age (in years) age+nhstepm*hstepm*stepm/12) by multiplying + nhstepm*hstepm matrices. + Output is stored in matrix po[i][j][h] for h every 'hstepm' step + (typically every 2 years instead of every month which is too big + for the memory). + Model is determined by parameters x and covariates have to be + included manually here. + + */ + + int i, j, d, h, k; + double **out, cov[NCOVMAX+1]; + double **newm; + double agexact; + double agebegin, ageend; + + /* Hstepm could be zero and should return the unit matrix */ + for (i=1;i<=nlstate+ndeath;i++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[i][j]=(i==j ? 1.0 : 0.0); + po[i][j][0]=(i==j ? 1.0 : 0.0); + } + /* Even if hstepm = 1, at least one multiplication by the unit matrix */ + for(h=1; h <=nhstepm; h++){ + for(d=1; d <=hstepm; d++){ + newm=savm; + /* Covariates have to be included here again */ + cov[1]=1.; + agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ + /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ + cov[2]=agexact; + if(nagesqr==1) + cov[3]= agexact*agexact; + for (k=1; k<=cptcovn;k++) + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ + /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ + for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ + + + /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ + /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ + out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, + oldm); + /* if((int)age == 70){ */ + /* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */ + /* for(i=1; i<=nlstate+ndeath; i++) { */ + /* printf("%d pmmij ",i); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",pmmij[i][j]); */ + /* } */ + /* printf(" oldm "); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",oldm[i][j]); */ + /* } */ + /* printf("\n"); */ + /* } */ + /* } */ + savm=oldm; + oldm=newm; + } + for(i=1; i<=nlstate+ndeath; i++) + for(j=1;j<=nlstate+ndeath;j++) { + po[i][j][h]=newm[i][j]; + /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/ + } + /*printf("h=%d ",h);*/ + } /* end h */ +/* printf("\n H=%d \n",h); */ + return po; +} + + +#ifdef NLOPT + double myfunc(unsigned n, const double *p1, double *grad, void *pd){ + double fret; + double *xt; + int j; + myfunc_data *d2 = (myfunc_data *) pd; +/* xt = (p1-1); */ + xt=vector(1,n); + for (j=1;j<=n;j++) xt[j]=p1[j-1]; /* xt[1]=p1[0] */ + + fret=(d2->function)(xt); /* p xt[1]@8 is fine */ + /* fret=(*func)(xt); /\* p xt[1]@8 is fine *\/ */ + printf("Function = %.12lf ",fret); + for (j=1;j<=n;j++) printf(" %d %.8lf", j, xt[j]); + printf("\n"); + free_vector(xt,1,n); + return fret; +} +#endif + +/*************** log-likelihood *************/ +double func( double *x) +{ + int i, ii, j, k, mi, d, kk; + double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; + double **out; + double sw; /* Sum of weights */ + double lli; /* Individual log likelihood */ + int s1, s2; + double bbh, survp; + long ipmx; double agexact; /*extern weight */ /* We are differentiating ll according to initial status */ @@ -2305,27 +2764,24 @@ double func( double *x) /* else */ /* lli=log(out[s1][s2] - savm[s1][s2]); */ /* #endif */ - lli=log(out[s1][s2] - savm[s1][s2]); - - } else if (s2==-2) { + lli=log(out[s1][s2] - savm[s1][s2]); + + } else if ( s2==-1 ) { /* alive */ 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 */ @@ -2437,6 +2893,10 @@ double func( double *x) s2=s[mw[mi+1][i]][i]; if( s2 > nlstate){ lli=log(out[s1][s2] - savm[s1][s2]); + } else if ( s2==-1 ) { /* alive */ + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= log(survp); }else{ lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ } @@ -2499,6 +2959,7 @@ double funcone( double *x) int s1, s2; double bbh, survp; double agexact; + double agebegin, ageend; /*extern weight */ /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ @@ -2517,7 +2978,12 @@ double funcone( double *x) oldm[ii][j]=(ii==j ? 1.0 : 0.0); savm[ii][j]=(ii==j ? 1.0 : 0.0); } - for(d=0; d nlstate && (mle <5) ){ /* Jackson */ lli=log(out[s1][s2] - savm[s1][s2]); - } else if (s2==-2) { + } else if ( s2==-1 ) { /* alive */ for (j=1,survp=0. ; j<=nlstate; j++) survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; lli= log(survp); @@ -2565,9 +3035,9 @@ double funcone( double *x) ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ if(globpr){ - fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f\ + fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ - num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], + num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); for(k=1,llt=0.,l=0.; k<=nlstate; k++){ llt +=ll[k]*gipmx/gsw; @@ -2605,8 +3075,8 @@ void likelione(FILE *ficres,double p[], printf("Problem with resultfile: %s\n", fileresilk); fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); } - fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); - fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight 2wlli out sav "); + fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); + fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav "); /* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ for(k=1; k<=nlstate; k++) fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); @@ -2616,11 +3086,23 @@ void likelione(FILE *ficres,double p[], *fretone=(*funcone)(p); if(*globpri !=0){ fclose(ficresilk); - fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle >= 1. You should at least run with mle >= 1 and starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); - fprintf(fichtm,"
- The first 3 individuals are drawn with lines. The function drawn is -2Log(L) in log scale: %s.png
\ -",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + if (mle ==0) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle = %d.",mle); + else if(mle >=1) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); + fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); + + + for (k=1; k<= nlstate ; k++) { + fprintf(fichtm,"
- Probability p%dj by origin %d and destination j. Dot's sizes are related to corresponding weight: %s-p%dj.png
\ +",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); + } + fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fprintf(fichtm,"
- and by state of destination %s-dest.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); fflush(fichtm); - } + } return; } @@ -2694,32 +3176,32 @@ void mlikeli(FILE *ficres,double p[], in #endif free_matrix(xi,1,npar,1,npar); fclose(ficrespow); - printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); - fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); } /**** Computes Hessian and covariance matrix ***/ -void hesscov(double **matcov, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) +void hesscov(double **matcov, double **hess, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) { double **a,**y,*x,pd; - double **hess; + /* double **hess; */ int i, j; int *indx; 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); + double hessij(double p[], double **hess, 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); + /* 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); - fprintf(ficlog,"%d",i);fflush(ficlog); + printf("%d-",i);fflush(stdout); + fprintf(ficlog,"%d-",i);fflush(ficlog); hess[i][i]=hessii(p,ftolhess,i,delti,func,npar); @@ -2730,9 +3212,9 @@ void hesscov(double **matcov, double p[] for (i=1;i<=npar;i++) { for (j=1;j<=npar;j++) { if (j>i) { - printf(".%d%d",i,j);fflush(stdout); - fprintf(ficlog,".%d%d",i,j);fflush(ficlog); - hess[i][j]=hessij(p,delti,i,j,func,npar); + printf(".%d-%d",i,j);fflush(stdout); + fprintf(ficlog,".%d-%d",i,j);fflush(ficlog); + hess[i][j]=hessij(p,hess, delti,i,j,func,npar); hess[j][i]=hess[i][j]; /*printf(" %lf ",hess[i][j]);*/ @@ -2766,53 +3248,78 @@ void hesscov(double **matcov, double p[] 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("%.6e ",hess[i][j]); + fprintf(ficlog,"%.6e ",hess[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } + /* printf("\n#Covariance matrix#\n"); */ + /* fprintf(ficlog,"\n#Covariance matrix#\n"); */ + /* for (i=1;i<=npar;i++) { */ + /* for (j=1;j<=npar;j++) { */ + /* printf("%.6e ",matcov[i][j]); */ + /* fprintf(ficlog,"%.6e ",matcov[i][j]); */ + /* } */ + /* printf("\n"); */ + /* fprintf(ficlog,"\n"); */ + /* } */ + /* Recompute Inverse */ - for (i=1;i<=npar;i++) - for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; - ludcmp(a,npar,indx,&pd); + /* for (i=1;i<=npar;i++) */ + /* for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; */ + /* ludcmp(a,npar,indx,&pd); */ + + /* printf("\n#Hessian matrix recomputed#\n"); */ + + /* for (j=1;j<=npar;j++) { */ + /* for (i=1;i<=npar;i++) x[i]=0; */ + /* x[j]=1; */ + /* lubksb(a,npar,indx,x); */ + /* 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"); */ + /* } */ + + /* Verifying the inverse matrix */ +#ifdef DEBUGHESS + y=matprod2(y,hess,1,npar,1,npar,1,npar,matcov); - /* printf("\n#Hessian matrix recomputed#\n"); + printf("\n#Verification: multiplying the matrix of covariance by the Hessian matrix, should be unity:#\n"); + fprintf(ficlog,"\n#Verification: multiplying the matrix of covariance by the Hessian matrix. Should be unity:#\n"); for (j=1;j<=npar;j++) { - for (i=1;i<=npar;i++) x[i]=0; - x[j]=1; - lubksb(a,npar,indx,x); for (i=1;i<=npar;i++){ - y[i][j]=x[i]; - printf("%.3e ",y[i][j]); - fprintf(ficlog,"%.3e ",y[i][j]); + printf("%.2f ",y[i][j]); + fprintf(ficlog,"%.2f ",y[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } - */ +#endif free_matrix(a,1,npar,1,npar); free_matrix(y,1,npar,1,npar); free_vector(x,1,npar); free_ivector(indx,1,npar); - free_matrix(hess,1,npar,1,npar); + /* free_matrix(hess,1,npar,1,npar); */ } /*************** hessian matrix ****************/ double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar) -{ +{ /* Around values of x, computes the function func and returns the scales delti and hessian */ int i; int l=1, lmax=20; - double k1,k2; + double k1,k2, res, fx; double p2[MAXPARM+1]; /* identical to x */ - double res; double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; - double fx; int k=0,kmax=10; double l1; @@ -2828,9 +3335,9 @@ double hessii(double x[], double delta, p2[theta]=x[theta]-delt; k2=func(p2)-fx; /*res= (k1-2.0*fx+k2)/delt/delt; */ - res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ + res= (k1+k2)/delt/delt/2.; /* Divided by 2 because L and not 2*L */ -#ifdef DEBUGHESS +#ifdef DEBUGHESSII 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 @@ -2844,48 +3351,125 @@ double hessii(double x[], double delta, else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ delts=delt; } - } + } /* End loop k */ } delti[theta]=delts; return res; } -double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) +double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,double (*func)(double []),int npar) { int i; int l=1, lmax=20; double k1,k2,k3,k4,res,fx; double p2[MAXPARM+1]; - int k; + int k, kmax=1; + double v1, v2, cv12, lc1, lc2; + int firstime=0; + fx=func(x); - for (k=1; k<=2; k++) { + for (k=1; k<=kmax; k=k+10) { for (i=1;i<=npar;i++) p2[i]=x[i]; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k1=func(p2)-fx; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k2=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k3=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k4=func(p2)-fx; - res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */ -#ifdef DEBUG - printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); - fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */ + if(k1*k2*k3*k4 <0.){ + firstime=1; + kmax=kmax+10; + } + if(kmax >=10 || firstime ==1){ + printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); + fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); + printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + } +#ifdef DEBUGHESSIJ + v1=hess[thetai][thetai]; + v2=hess[thetaj][thetaj]; + cv12=res; + /* Computing eigen value of Hessian matrix */ + lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + if ((lc2 <0) || (lc1 <0) ){ + printf("Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + } #endif } return res; } + /* Not done yet: Was supposed to fix if not exactly at the maximum */ +/* double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) */ +/* { */ +/* int i; */ +/* int l=1, lmax=20; */ +/* double k1,k2,k3,k4,res,fx; */ +/* double p2[MAXPARM+1]; */ +/* double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; */ +/* int k=0,kmax=10; */ +/* double l1; */ + +/* fx=func(x); */ +/* for(l=0 ; l <=lmax; l++){ /\* Enlarging the zone around the Maximum *\/ */ +/* l1=pow(10,l); */ +/* delts=delt; */ +/* for(k=1 ; k khi/nkhif) || (k2 >khi/nkhif) || (k4 >khi/nkhif) || (k4 >khi/nkhif)){ /\* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. *\/ */ +/* k=kmax; l=lmax*10; */ +/* } */ +/* else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ */ +/* delts=delt; */ +/* } */ +/* } /\* End loop k *\/ */ +/* } */ +/* delti[theta]=delts; */ +/* return res; */ +/* } */ + + /************** Inverse of matrix **************/ void ludcmp(double **a, int n, int *indx, double *d) { @@ -2966,25 +3550,61 @@ void pstamp(FILE *fichier) } /************ Frequencies ********************/ -void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) +void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ + int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],\ + int firstpass, int lastpass, int stepm, int weightopt, char model[]) { /* Some frequencies */ int i, m, jk, j1, bool, z1,j; + int mi; /* Effective wave */ int first; double ***freq; /* Frequencies */ double *pp, **prop; double pos,posprop, k2, dateintsum=0,k2cpt=0; - char fileresp[FILENAMELENGTH]; - + char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; + double agebegin, ageend; + pp=vector(1,nlstate); prop=matrix(1,nlstate,iagemin,iagemax+3); strcpy(fileresp,"P_"); strcat(fileresp,fileresu); + /*strcat(fileresphtm,fileresu);*/ if((ficresp=fopen(fileresp,"w"))==NULL) { printf("Problem with prevalence resultfile: %s\n", fileresp); fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } + + strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); + if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { + printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtm,"\nIMaCh PHTM_ %s\n %s
%s \ +
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ + fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies and prevalence by age at begin of transition

\n",fileresphtm, fileresphtm); + + strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); + if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { + printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s \ +
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ + fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions by age at begin of transition

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); j1=0; @@ -2993,10 +3613,7 @@ void freqsummary(char fileres[], int ia first=1; - /* for(k1=1; k1<=j ; k1++){ */ /* Loop on covariates */ - /* for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */ - /* j1++; */ - for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ + for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */ /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ for (i=-5; i<=nlstate+ndeath; i++) @@ -3010,9 +3627,9 @@ void freqsummary(char fileres[], int ia dateintsum=0; k2cpt=0; - for (i=1; i<=imx; i++) { + for (i=1; i<=imx; i++) { /* For each individual i */ bool=1; - if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ for (z1=1; z1<=cptcoveff; z1++) if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* Tests if the value of each of the covariates of i is equal to filter j1 */ @@ -3022,51 +3639,98 @@ void freqsummary(char fileres[], int ia j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ } - } - + } /* cptcovn > 0 */ + if (bool==1){ - for(m=firstpass; m<=lastpass; m++){ - k2=anint[m][i]+(mint[m][i]/12.); - /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ - if(agev[m][i]==0) agev[m][i]=iagemax+1; - if(agev[m][i]==1) agev[m][i]=iagemax+2; - if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; + /* for(m=firstpass; m<=lastpass; m++){ */ + for(mi=1; mi=firstpass && m <=lastpass){ + k2=anint[m][i]+(mint[m][i]/12.); + /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ + if(agev[m][i]==0) agev[m][i]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ + if(agev[m][i]==1) agev[m][i]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ + if (s[m][i]>0 && s[m][i]<=nlstate) /* If status at wave m is known and a live state */ + prop[s[m][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */ if (m1) && (agev[m][i]< (iagemax+3))) { - dateintsum=dateintsum+k2; - k2cpt++; - } - /*}*/ - } - } - } /* end i */ + } + if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) { + dateintsum=dateintsum+k2; + k2cpt++; + /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ + } + /*}*/ + } /* end m */ + } /* end bool */ + } /* end i = 1 to imx */ /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ pstamp(ficresp); if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "\n

********** Variable "); + fprintf(ficresphtmfr, "\n

********** Variable "); + for (z1=1; z1<=cptcoveff; z1++){ + fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + } + fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "**********

\n"); + fprintf(ficresphtmfr, "**********\n"); fprintf(ficlog, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficlog, "**********\n#"); + fprintf(ficlog, "**********\n"); } - for(i=1; i<=nlstate;i++) + fprintf(ficresphtm,""); + for(i=1; i<=nlstate;i++) { fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); + fprintf(ficresphtm, "",i,i); + } fprintf(ficresp, "\n"); + fprintf(ficresphtm, "\n"); + /* Header of frequency table by age */ + fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); + fprintf(ficresphtmfr," "); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",jk,m); + } + } + fprintf(ficresphtmfr, "\n"); + + /* For each age */ for(i=iagemin; i <= iagemax+3; i++){ - if(i==iagemax+3){ + fprintf(ficresphtm,""); + if(i==iagemax+1){ + fprintf(ficlog,"1"); + fprintf(ficresphtmfr," "); + }else if(i==iagemax+2){ + fprintf(ficlog,"0"); + fprintf(ficresphtmfr," "); + }else if(i==iagemax+3){ fprintf(ficlog,"Total"); + fprintf(ficresphtmfr," "); }else{ if(first==1){ first=0; printf("See log file for details...\n"); } + fprintf(ficresphtmfr," ",i); fprintf(ficlog,"Age %d", i); } for(jk=1; jk <=nlstate ; jk++){ @@ -3109,32 +3773,47 @@ void freqsummary(char fileres[], int ia if( i <= iagemax){ if(pos>=1.e-5){ fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); + fprintf(ficresphtm,"",i,prop[jk][i]/posprop, prop[jk][i],posprop); /*probs[i][jk][j1]= pp[jk]/pos;*/ /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ } - else + else{ fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); + fprintf(ficresphtm,"",i, prop[jk][i],posprop); + } } } - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) - if(freq[jk][m][i] !=0 ) { - if(first==1) - printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[jk][m][i] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); + } fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]); } - if(i <= iagemax) + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",freq[jk][m][i]); + } + } + fprintf(ficresphtmfr,"\n "); + if(i <= iagemax){ fprintf(ficresp,"\n"); + fprintf(ficresphtm,"\n"); + } if(first==1) printf("Others in log...\n"); fprintf(ficlog,"\n"); - } + } /* end loop i */ + fprintf(ficresphtm,"
Age%d%d
0
Unknown
Total
%d%d%.5f%.0f%.0f%dNaNq%.0f%.0f%.0f
\n"); + fprintf(ficresphtmfr,"\n"); /*}*/ - } + } /* end j1 */ dateintmean=dateintsum/k2cpt; fclose(ficresp); + fclose(ficresphtm); + fclose(ficresphtmfr); free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); free_vector(pp,1,nlstate); free_matrix(prop,1,nlstate,iagemin, iagemax+3); @@ -3150,6 +3829,9 @@ void prevalence(double ***probs, double */ int i, m, jk, j1, bool, z1,j; + int mi; /* Effective wave */ + int iage; + double agebegin, ageend; double **prop; double posprop; @@ -3169,54 +3851,57 @@ void prevalence(double ***probs, double first=1; for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ - /*for(i1=1; i1<=ncodemax[k1];i1++){ - j1++;*/ - - for (i=1; i<=nlstate; i++) - for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0.0; - - for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) - bool=0; - } - if (bool==1) { - for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/ + for (i=1; i<=nlstate; i++) + for(iage=iagemin; iage <= iagemax+3; iage++) + prop[i][iage]=0.0; + + for (i=1; i<=imx; i++) { /* Each individual */ + bool=1; + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcoveff; z1++) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) + bool=0; + } + if (bool==1) { + /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ + for(mi=1; mi=firstpass && m <=lastpass){ y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ if(agev[m][i]==0) agev[m][i]=iagemax+1; if(agev[m][i]==1) agev[m][i]=iagemax+2; if((int)agev[m][i] iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); - if (s[m][i]>0 && s[m][i]<=nlstate) { + if (s[m][i]>0 && s[m][i]<=nlstate) { /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ - prop[s[m][i]][(int)agev[m][i]] += weight[i]; - prop[s[m][i]][iagemax+3] += weight[i]; - } - } + prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ + prop[s[m][i]][iagemax+3] += weight[i]; + } /* end valid statuses */ + } /* end selection of dates */ } /* end selection of waves */ - } - } - for(i=iagemin; i <= iagemax+3; i++){ - for(jk=1,posprop=0; jk <=nlstate ; jk++) { - posprop += prop[jk][i]; - } - - for(jk=1; jk <=nlstate ; jk++){ - if( i <= iagemax){ - if(posprop>=1.e-5){ - probs[i][jk][j1]= prop[jk][i]/posprop; - } else{ - if(first==1){ - first=0; - printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); - } + } /* end effective waves */ + } /* end bool */ + } + for(i=iagemin; i <= iagemax+3; i++){ + for(jk=1,posprop=0; jk <=nlstate ; jk++) { + posprop += prop[jk][i]; + } + + for(jk=1; jk <=nlstate ; jk++){ + if( i <= iagemax){ + if(posprop>=1.e-5){ + probs[i][jk][j1]= prop[jk][i]/posprop; + } else{ + if(first==1){ + first=0; + printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); } - } - }/* end jk */ - }/* end i */ + } + } + }/* end jk */ + }/* end i */ /*} *//* end i1 */ } /* end j1 */ @@ -3239,31 +3924,67 @@ void concatwav(int wav[], int **dh, int int i, mi, m; /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; double sum=0., jmean=0.;*/ - int first; + int first, firstwo, firsthree; int j, k=0,jk, ju, jl; double sum=0.; first=0; + firstwo=0; + firsthree=0; jmin=100000; jmax=-1; jmean=0.; - for(i=1; i<=imx; i++){ + for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ mi=0; m=firstpass; - while(s[m][i] <= nlstate){ - if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5) + while(s[m][i] <= nlstate){ /* a live state */ + if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ mw[++mi][i]=m; - if(m >=lastpass) + } + if(m >=lastpass){ + if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ + if(firsthree == 0){ + printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + firsthree=1; + }else{ + fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + } + mw[++mi][i]=m; + } + if(s[m][i]==-2){ /* Vital status is really unknown */ + nbwarn++; + if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ + printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); + fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); + } + break; + } break; + } else m++; }/* end while */ - if (s[m][i] > nlstate){ + + /* After last pass */ + if (s[m][i] > nlstate){ /* In a death state */ mi++; /* Death is another wave */ /* if(mi==0) never been interviewed correctly before death */ /* Only death is a correct wave */ mw[mi][i]=m; + }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */ + /* m++; */ + /* mi++; */ + /* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */ + /* mw[mi][i]=m; */ + nberr++; + if(firstwo==0){ + printf("Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + firstwo=1; + }else if(firstwo==1){ + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + } } - wav[i]=mi; if(mi==0){ nbwarn++; @@ -3276,7 +3997,9 @@ void concatwav(int wav[], int **dh, int } } /* end mi==0 */ } /* End individuals */ + /* wav and mw are no more changed */ + for(i=1; i<=imx; i++){ for(mi=1; mi= YEARM) hstepm=1; nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ gradg=matrix(1,npar,1,nlstate); + mgp=matrix(1,npar,1,nlstate); + mgm=matrix(1,npar,1,nlstate); gp=vector(1,nlstate); gm=vector(1,nlstate); @@ -4167,18 +4897,27 @@ void varprevlim(char fileres[], double * for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - for(i=1;i<=nlstate;i++) + if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + else + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + for(i=1;i<=nlstate;i++){ gp[i] = prlim[i][i]; - + mgp[theta][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,ij); - for(i=1;i<=nlstate;i++) + if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + else + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + for(i=1;i<=nlstate;i++){ gm[i] = prlim[i][i]; - + mgm[theta][i] = prlim[i][i]; + } for(i=1;i<=nlstate;i++) gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; + /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */ } /* End theta */ trgradg =matrix(1,nlstate,1,npar); @@ -4186,11 +4925,34 @@ void varprevlim(char fileres[], double * for(j=1; j<=nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[j][theta]=gradg[theta][j]; + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\nmgm mgp %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf(" %d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\n gradg %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf("%d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf("%d %lf ",theta,gradg[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ for(i=1;i<=nlstate;i++) varpl[i][(int)age] =0.; + if((int)age==79 ||(int)age== 80 ||(int)age== 81){ + matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + }else{ matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + } for(i=1;i<=nlstate;i++) varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ @@ -4200,6 +4962,8 @@ void varprevlim(char fileres[], double * fprintf(ficresvpl,"\n"); free_vector(gp,1,nlstate); free_vector(gm,1,nlstate); + free_matrix(mgm,1,npar,1,nlstate); + free_matrix(mgp,1,npar,1,nlstate); free_matrix(gradg,1,npar,1,nlstate); free_matrix(trgradg,1,nlstate,1,npar); } /* End age */ @@ -4550,30 +5314,41 @@ To be simple, these graphs help to under void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ int lastpass, int stepm, int weightopt, char model[],\ int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ - int popforecast, int estepm ,\ - double jprev1, double mprev1,double anprev1, \ - double jprev2, double mprev2,double anprev2){ + int popforecast, int prevfcast, int backcast, int estepm , \ + double jprev1, double mprev1,double anprev1, double dateprev1, \ + double jprev2, double mprev2,double anprev2, double dateprev2){ int jj1, k1, i1, cpt; fprintf(fichtm,""); - fprintf(fichtm,"