--- imach/src/imach.c 2015/04/23 12:01:52 1.186 +++ imach/src/imach.c 2015/09/07 14:09:23 1.199 @@ -1,6 +1,52 @@ -/* $Id: imach.c,v 1.186 2015/04/23 12:01:52 brouard Exp $ +/* $Id: imach.c,v 1.199 2015/09/07 14:09:23 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.199 2015/09/07 14:09:23 brouard + Summary: 0.98q6 changing default small png format for graph to vectorized svg. + + Revision 1.198 2015/09/03 07:14:39 brouard + Summary: 0.98q5 Flavia + + Revision 1.197 2015/09/01 18:24:39 brouard + *** empty log message *** + + Revision 1.196 2015/08/18 23:17:52 brouard + Summary: 0.98q5 + + Revision 1.195 2015/08/18 16:28:39 brouard + Summary: Adding a hack for testing purpose + + After reading the title, ftol and model lines, if the comment line has + a q, starting with #q, the answer at the end of the run is quit. It + permits to run test files in batch with ctest. The former workaround was + $ echo q | imach foo.imach + + Revision 1.194 2015/08/18 13:32:00 brouard + Summary: Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line. + + Revision 1.193 2015/08/04 07:17:42 brouard + Summary: 0.98q4 + + Revision 1.192 2015/07/16 16:49:02 brouard + Summary: Fixing some outputs + + Revision 1.191 2015/07/14 10:00:33 brouard + Summary: Some fixes + + Revision 1.190 2015/05/05 08:51:13 brouard + Summary: Adding digits in output parameters (7 digits instead of 6) + + Fix 1+age+. + + Revision 1.189 2015/04/30 14:45:16 brouard + Summary: 0.98q2 + + Revision 1.188 2015/04/30 08:27:53 brouard + *** empty log message *** + + Revision 1.187 2015/04/29 09:11:15 brouard + *** empty log message *** + Revision 1.186 2015/04/23 12:01:52 brouard Summary: V1*age is working now, version 0.98q1 @@ -585,7 +631,10 @@ end */ +/* #define DEBUG */ +/* #define DEBUGBRENT */ #define POWELL /* Instead of NLOPT */ +#define POWELLF1F3 /* Skip test */ /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ @@ -654,11 +703,12 @@ typedef struct { #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ #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) ; +#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 #define MAXN 20000 #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 #define AGEBASE 40 +#define AGEOVERFLOW 1.e20 #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */ #ifdef _WIN32 #define DIRSEPARATOR '\\' @@ -670,15 +720,16 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.186 2015/04/23 12:01:52 brouard Exp $ */ +/* $Id: imach.c,v 1.199 2015/09/07 14:09:23 brouard Exp $ */ /* $State: Exp $ */ - -char version[]="Imach version 0.98q1, April 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.186 $ $Date: 2015/04/23 12:01:52 $"; +#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.199 $ $Date: 2015/09/07 14:09:23 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ -int nvar=0, nforce=0; /* Number of variables, number of forces */ +int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */ /* Number of covariates model=V2+V1+ V3*age+V2*V4 */ int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */ int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */ @@ -808,7 +859,13 @@ int estepm; int m,nb; long *num; -int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage,*cens; +int firstpass=0, lastpass=4,*cod, *cens; +int *ncodemax; /* ncodemax[j]= Number of modalities of the j th + covariate for which somebody answered excluding + undefined. Usually 2: 0 and 1. */ +int *ncodemaxwundef; /* ncodemax[j]= Number of modalities of the j th + covariate for which somebody answered including + undefined. Usually 3: -1, 0 and 1. */ double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; double **pmmij, ***probs; double *ageexmed,*agecens; @@ -819,9 +876,10 @@ int **s; /* Status */ double *agedc; double **covar; /**< covar[j,i], value of jth covariate for individual i, * covar=matrix(0,NCOVMAX,1,n); - * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; */ + * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ double idx; int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ +int *Tage; int *Ndum; /** Freq of modality (tricode */ int **codtab; /**< codtab=imatrix(1,100,1,10); */ int **Tvard, *Tprod, cptcovprod, *Tvaraff; @@ -915,11 +973,54 @@ char *trimbb(char *out, char *in) return s; } +/* char *substrchaine(char *out, char *in, char *chain) */ +/* { */ +/* /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */ +/* char *s, *t; */ +/* t=in;s=out; */ +/* while ((*in != *chain) && (*in != '\0')){ */ +/* *out++ = *in++; */ +/* } */ + +/* /\* *in matches *chain *\/ */ +/* while ((*in++ == *chain++) && (*in != '\0')){ */ +/* printf("*in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ +/* } */ +/* in--; chain--; */ +/* while ( (*in != '\0')){ */ +/* printf("Bef *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ +/* *out++ = *in++; */ +/* printf("Aft *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ +/* } */ +/* *out='\0'; */ +/* out=s; */ +/* return out; */ +/* } */ +char *substrchaine(char *out, char *in, char *chain) +{ + /* Substract chain 'chain' from 'in', return and output 'out' */ + /* in="V1+V1*age+age*age+V2", chain="age*age" */ + + char *strloc; + + strcpy (out, in); + strloc = strstr(out, chain); /* strloc points to out at age*age+V2 */ + printf("Bef strloc=%s chain=%s out=%s \n", strloc, chain, out); + if(strloc != NULL){ + /* will affect out */ /* strloc+strlenc(chain)=+V2 */ /* Will also work in Unicode */ + memmove(strloc,strloc+strlen(chain), strlen(strloc+strlen(chain))+1); + /* strcpy (strloc, strloc +strlen(chain));*/ + } + printf("Aft strloc=%s chain=%s in=%s out=%s \n", strloc, chain, in, out); + return out; +} + + char *cutl(char *blocc, char *alocc, char *in, char occ) { - /* cuts string in into blocc and alocc where blocc ends before first occurence of char 'occ' + /* cuts string in into blocc and alocc where blocc ends before FIRST occurence of char 'occ' and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') - gives blocc="abcdef2ghi" and alocc="j". + gives blocc="abcdef" and alocc="ghi2j". If occ is not found blocc is null and alocc is equal to in. Returns blocc */ char *s, *t; @@ -945,7 +1046,7 @@ char *cutl(char *blocc, char *alocc, cha } char *cutv(char *blocc, char *alocc, char *in, char occ) { - /* cuts string in into blocc and alocc where blocc ends before last occurence of char 'occ' + /* cuts string in into blocc and alocc where blocc ends before LAST occurence of char 'occ' and alocc starts after last occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') gives blocc="abcdef2ghi" and alocc="j". If occ is not found blocc is null and alocc is equal to in. Returns alocc @@ -1256,7 +1357,13 @@ double f1dim(double x) /*****************brent *************************/ double brent(double ax, double bx, double cx, double (*f)(double), double tol, double *xmin) -{ +{ + /* Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is + * between ax and cx, and f(bx) is less than both f(ax) and f(cx) ), this routine isolates + * the minimum to a fractional precision of about tol using Brent’s method. The abscissa of + * the minimum is returned as xmin, and the minimum function value is returned as brent , the + * returned function value. + */ int iter; double a,b,d,etemp; double fu=0,fv,fw,fx; @@ -1339,9 +1446,20 @@ values at the three points, fa, fb , and */ double ulim,u,r,q, dum; double fu; - - *fa=(*func)(*ax); - *fb=(*func)(*bx); + + double scale=10.; + int iterscale=0; + + *fa=(*func)(*ax); /* xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/ + *fb=(*func)(*bx); /* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */ + + + /* while(*fb != *fb){ /\* *ax should be ok, reducing distance to *ax *\/ */ + /* printf("Warning mnbrak *fb = %lf, *bx=%lf *ax=%lf *fa==%lf iter=%d\n",*fb, *bx, *ax, *fa, iterscale++); */ + /* *bx = *ax - (*ax - *bx)/scale; */ + /* *fb=(*func)(*bx); /\* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) *\/ */ + /* } */ + if (*fb > *fa) { SHFT(dum,*ax,*bx,dum) SHFT(dum,*fb,*fa,dum) @@ -1374,31 +1492,41 @@ values at the three points, fa, fb , and #endif #ifdef MNBRAKORIGINAL #else - if (fu > *fc) { -#ifdef DEBUG - printf("mnbrak4 fu > fc \n"); - fprintf(ficlog, "mnbrak4 fu > fc\n"); -#endif - /* SHFT(u,*cx,*cx,u) /\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\/ */ - /* SHFT(*fa,*fc,fu,*fc) /\* (b, u, c) is a bracket while test fb > fc will be fu > fc will exit *\/ */ - dum=u; /* Shifting c and u */ - u = *cx; - *cx = dum; - dum = fu; - fu = *fc; - *fc =dum; - } else { /* end */ +/* if (fu > *fc) { */ +/* #ifdef DEBUG */ +/* printf("mnbrak4 fu > fc \n"); */ +/* fprintf(ficlog, "mnbrak4 fu > fc\n"); */ +/* #endif */ +/* /\* SHFT(u,*cx,*cx,u) /\\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\\/ *\/ */ +/* /\* SHFT(*fa,*fc,fu,*fc) /\\* (b, u, c) is a bracket while test fb > fc will be fu > fc will exit *\\/ *\/ */ +/* dum=u; /\* Shifting c and u *\/ */ +/* u = *cx; */ +/* *cx = dum; */ +/* dum = fu; */ +/* fu = *fc; */ +/* *fc =dum; */ +/* } else { /\* end *\/ */ +/* #ifdef DEBUG */ +/* printf("mnbrak3 fu < fc \n"); */ +/* fprintf(ficlog, "mnbrak3 fu < fc\n"); */ +/* #endif */ +/* dum=u; /\* Shifting c and u *\/ */ +/* u = *cx; */ +/* *cx = dum; */ +/* dum = fu; */ +/* fu = *fc; */ +/* *fc =dum; */ +/* } */ #ifdef DEBUG - printf("mnbrak3 fu < fc \n"); - fprintf(ficlog, "mnbrak3 fu < fc\n"); + printf("mnbrak34 fu < or >= fc \n"); + fprintf(ficlog, "mnbrak34 fu < fc\n"); #endif - dum=u; /* Shifting c and u */ - u = *cx; - *cx = dum; - dum = fu; - fu = *fc; - *fc =dum; - } + dum=u; /* Shifting c and u */ + u = *cx; + *cx = dum; + dum = fu; + fu = *fc; + *fc =dum; #endif } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */ #ifdef DEBUG @@ -1458,6 +1586,8 @@ void linmin(double p[], double xi[], int int j; double xx,xmin,bx,ax; double fx,fb,fa; + + double scale=10., axs, xxs, xxss; /* Scale added for infinity */ ncom=n; pcom=vector(1,n); @@ -1467,18 +1597,57 @@ void linmin(double p[], double xi[], int pcom[j]=p[j]; xicom[j]=xi[j]; } - ax=0.0; - xx=1.0; - mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */ - *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */ + + /* axs=0.0; */ + /* xxss=1; /\* 1 and using scale *\/ */ + xxs=1; + /* do{ */ + ax=0.; + xx= xxs; + 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)) */ + /* Outputs: fa=f(p(j)) and fx=f(p(j) + xxs * xi(j) ) and f(bx)= f(p(j)+ bx* xi(j)) */ + /* 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 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); +#endif + *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/ + /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */ + /* fmin = f(p[j] + xmin * xi[j]) */ + /* P+lambda n in that direction (lambdamin), with TOL between abscisses */ + /* f1dim(xmin): for (j=1;j<=ncom;j++) xt[j]=pcom[j]+xmin*xicom[j]; */ #ifdef DEBUG printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); #endif +#ifdef DEBUGLINMIN + printf("linmin end "); +#endif for (j=1;j<=n;j++) { - xi[j] *= xmin; - p[j] += xi[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 ); */ + p[j] += xi[j]; /* Parameters values are updated accordingly */ } + /* printf("\n"); */ +#ifdef DEBUGLINMIN + printf("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++) { + printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]); + if(j % ncovmodel == 0) + printf("\n"); + } +#endif free_vector(xicom,1,n); free_vector(pcom,1,n); } @@ -1513,7 +1682,7 @@ void powell(double p[], double **xi, int for (j=1;j<=n;j++) pt[j]=p[j]; rcurr_time = time(NULL); for (*iter=1;;++(*iter)) { - fp=(*fret); + fp=(*fret); /* From former iteration or initial value */ ibig=0; del=0.0; rlast_time=rcurr_time; @@ -1523,7 +1692,7 @@ void powell(double p[], double **xi, int printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); /* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ - for (i=1;i<=n;i++) { + for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); fprintf(ficrespow," %.12lf", p[i]); @@ -1551,22 +1720,22 @@ void powell(double p[], double **xi, int 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(rforecast_time-rcurr_time,tmpout),strfor,strcurr); } } - for (i=1;i<=n;i++) { - for (j=1;j<=n;j++) xit[j]=xi[j][i]; + for (i=1;i<=n;i++) { /* For each direction i */ + 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); #endif - printf("%d",i);fflush(stdout); + printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); - linmin(p,xit,n,fret,func); /* xit[n] has been loaded for direction i */ - if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions - because that direction will be replaced unless the gain del is small - in comparison with the 'probable' gain, mu^2, with the last average direction. - Unless the n directions are conjugate some gain in the determinant may be obtained - with the new direction. - */ + 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 */ + if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */ + /* because that direction will be replaced unless the gain del is small */ + /* in comparison with the 'probable' gain, mu^2, with the last average direction. */ + /* Unless the n directions are conjugate some gain in the determinant may be obtained */ + /* with the new direction. */ del=fabs(fptt-(*fret)); ibig=i; } @@ -1585,8 +1754,23 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } /* end i */ + } /* end loop on each direction i */ + /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ + /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ + /* New value of last point Pn is not computed, P(n-1) */ if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */ + /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ + /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */ + /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */ + /* decreased of more than 3.84 */ + /* By adding age*age and V1*age the gain (-2LL) should be more than 5.99 (ddl=2) */ + /* By using V1+V2+V3, the gain should be 7.82, compared with basic 1+age. */ + /* By adding 10 parameters more the gain should be 18.31 */ + + /* Starting the program with initial values given by a former maximization will simply change */ + /* the scales of the directions and the directions, because the are reset to canonical directions */ + /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */ + /* under the tolerance value. If the tolerance is very small 1.e-9, it could last long. */ #ifdef DEBUG int k[2],l; k[0]=1; @@ -1616,7 +1800,7 @@ void powell(double p[], double **xi, int free_vector(ptt,1,n); free_vector(pt,1,n); return; - } + } /* enough precision */ if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ ptt[j]=2.0*p[j]-pt[j]; @@ -1624,7 +1808,10 @@ void powell(double p[], double **xi, int pt[j]=p[j]; } fptt=(*func)(ptt); /* f_3 */ +#ifdef POWELLF1F3 +#else if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */ +#endif /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */ /* From x1 (P0) distance of x2 is at h and x3 is 2h */ /* Let f"(x2) be the 2nd derivative equal everywhere. */ @@ -1653,14 +1840,29 @@ void powell(double p[], double **xi, int if (t < 0.0) { /* Then we use it for new direction */ #else if (directest*t < 0.0) { /* Contradiction between both tests */ - printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); - printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); - fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); - fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); - } + printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); + printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); + fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); + fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); + } if (directest < 0.0) { /* Then we use it for new direction */ #endif - linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/ +#ifdef DEBUGLINMIN + printf("Before linmin in direction P%d-P0\n",n); + for (j=1;j<=n;j++) { + printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0) + printf("\n"); + } +#endif + linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ +#ifdef DEBUGLINMIN + for (j=1;j<=n;j++) { + printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0) + printf("\n"); + } +#endif for (j=1;j<=n;j++) { xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ @@ -1678,9 +1880,12 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } /* end of t negative */ + } /* end of t or directest negative */ +#ifdef POWELLF1F3 +#else } /* end if (fptt < fp) */ - } +#endif + } /* loop iteration */ } /**** Prevalence limit (stable or period prevalence) ****************/ @@ -1709,15 +1914,16 @@ double **prevalim(double **prlim, int nl 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+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; - /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/ + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[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+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[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<=cptcovprod;k++) /* Useless */ - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(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]);*/ @@ -1871,6 +2077,7 @@ double ***hpxij(double ***po, int nhstep int i, j, d, h, k; double **out, cov[NCOVMAX+1]; double **newm; + double agexact; /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) @@ -1884,14 +2091,17 @@ double ***hpxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM; + agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; + cov[2]=agexact; + if(nagesqr==1) + cov[3]= agexact*agexact; for (k=1; k<=cptcovn;k++) - cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[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+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[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+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + 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);*/ @@ -1943,6 +2153,7 @@ double func( double *x) int s1, s2; double bbh, survp; long ipmx; + double agexact; /*extern weight */ /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ @@ -1964,7 +2175,7 @@ double func( double *x) to be observed in j being in i according to the model. */ for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */ - cov[2+k]=covar[Tvar[k]][i]; + cov[2+nagesqr+k]=covar[Tvar[k]][i]; } /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] @@ -1977,9 +2188,12 @@ double func( double *x) } for(d=0; d0) { /* 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]][codtab[j1][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 */ bool=0; - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", - bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1], - j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/ - /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/ + /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", + bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), + 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*/ } } @@ -2776,10 +3007,10 @@ void freqsummary(char fileres[], int ia 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]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresp, "**********\n#"); fprintf(ficlog, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficlog, "**********\n#"); } for(i=1; i<=nlstate;i++) @@ -2907,7 +3138,7 @@ void prevalence(double ***probs, double bool=1; if (cptcovn>0) { for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) bool=0; } if (bool==1) { @@ -3112,11 +3343,11 @@ void tricode(int *Tvar, int **nbcode, in cptcoveff=0; - for (k=-1; k < maxncov; k++) Ndum[k]=0; for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ /* Loop on covariates without age and products */ for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */ + for (k=-1; k < maxncov; k++) Ndum[k]=0; for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i @@ -3139,15 +3370,24 @@ void tricode(int *Tvar, int **nbcode, in /* getting the maximum value of the modality of the covariate (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and female is 1, then modmaxcovj=1.*/ - } + } /* end for loop on individuals i */ printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); + fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); cptcode=modmaxcovj; /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ /*for (i=0; i<=cptcode; i++) {*/ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */ - printf("Frequencies of covariates %d V%d %d\n", j, Tvar[j], Ndum[i]); - if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */ - ncodemax[j]++; /* ncodemax[j]= Number of non-null modalities of the j th covariate. */ + for (k=modmincovj; k<=modmaxcovj; k++) { /* k=-1 ? 0 and 1*//* For each value k of the modality of model-cov j */ + printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); + fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); + if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */ + if( k != -1){ + ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th + covariate for which somebody answered excluding + undefined. Usually 2: 0 and 1. */ + } + ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th + covariate for which somebody answered including + undefined. Usually 3: -1, 0 and 1. */ } /* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ @@ -3164,41 +3404,53 @@ void tricode(int *Tvar, int **nbcode, in nbcode[Tvar[j]][1]=0; nbcode[Tvar[j]][2]=1; nbcode[Tvar[j]][3]=2; + To be continued (not working yet). */ - ij=1; /* ij is similar to i but can jumps over null modalities */ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 */ - for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */ - /*recode from 0 */ - if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */ - nbcode[Tvar[j]][ij]=k; /* stores the modality k in an array nbcode. - k is a modality. If we have model=V1+V1*sex - then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ - ij++; - } - if (ij > ncodemax[j]) break; - } /* end of loop on */ - } /* end of loop on modality */ + ij=0; /* ij is similar to i but can jump over null modalities */ + for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + if (Ndum[i] == 0) { /* If nobody responded to this modality k */ + break; + } + ij++; + nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ + cptcode = ij; /* New max modality for covar j */ + } /* end of loop on modality i=-1 to 1 or more */ + + /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ + /* /\*recode from 0 *\/ */ + /* k is a modality. If we have model=V1+V1*sex */ + /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ + /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ + /* } */ + /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ + /* if (ij > ncodemax[j]) { */ + /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ + /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ + /* break; */ + /* } */ + /* } /\* end of loop on modality k *\/ */ } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */ + for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ - Ndum[ij]++; + Ndum[ij]++; /* Might be supersed V1 + V1*age */ } - ij=1; + ij=0; for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ if((Ndum[i]!=0) && (i<=ncovcol)){ + ij++; /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ Tvaraff[ij]=i; /*For printing (unclear) */ - ij++; - }else - Tvaraff[ij]=0; + }else{ + /* Tvaraff[ij]=0; */ + } } - ij--; + /* ij--; */ cptcoveff=ij; /*Number of total covariates*/ } @@ -3802,7 +4054,8 @@ void varevsij(char optionfilefiname[], d 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 size 320, 240"); + /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */ + fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480"); /* 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); */ @@ -3812,11 +4065,11 @@ void varevsij(char optionfilefiname[], d fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev)); fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 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(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); +/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */ + fprintf(ficgp,"\nset out \"%s%s.svg\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,nlstate); @@ -3989,10 +4242,9 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Matrix of variance-covariance of pairs of step probabilities (drawings)

    • \n",optionfilehtmcov); - fprintf(fichtmcov,"\n

    Matrix of variance-covariance of pairs of step probabilities

    \n\ - file %s
    \n",optionfilehtmcov); - fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated\ + fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)


    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery
    • \n",optionfilehtmcov); + fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); + fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ and drawn. It helps understanding how is the covariance between two incidences.\ They are expressed in year-1 in order to be less dependent of stepm.
    \n"); fprintf(fichtmcov,"\n
    Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \ @@ -4013,23 +4265,23 @@ To be simple, these graphs help to under /*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]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(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]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(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]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(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]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(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]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresprobcor, "**********\n#"); } @@ -4039,8 +4291,10 @@ To be simple, these graphs help to under gm=vector(1,(nlstate)*(nlstate+ndeath)); for (age=bage; age<=fage; age ++){ cov[2]=age; + if(nagesqr==1) + cov[3]= age*age; for (k=1; k<=cptcovn;k++) { - cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4 + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];/* j1 1 2 3 4 * 1 1 1 1 1 * 2 2 1 1 1 * 3 1 2 1 1 @@ -4048,9 +4302,9 @@ To be simple, these graphs help to under /* nbcode[1][1]=0 nbcode[1][2]=1;*/ } /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; + for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[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]]]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; for(theta=1; theta <=npar; theta++){ @@ -4200,15 +4454,15 @@ To be simple, these graphs help to under 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 size 320, 240"); + fprintf(ficgp,"\nset ter svg size 640, 480"); 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,\ + :\ +%s%d%1d%1d-%1d%1d.svg, ",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
    ",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 out \"%s%d%1d%1d-%1d%1d.svg\"",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",\ @@ -4225,7 +4479,7 @@ To be simple, these graphs help to under }/* 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); + fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.svg\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); first=1; } /*l12 */ } /* k12 */ @@ -4282,40 +4536,49 @@ fprintf(fichtm," \n"); - fprintf(fichtm,"\ \n

  • Result files (second order: variances)

    \n\ - - Parameter file with estimated parameters and covariance matrix: %s
    \n", rfileres,rfileres); + - Parameter file with estimated parameters and covariance matrix: %s
    \ + - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.
    \ +But because parameters are usually highly correlated (a higher incidence of disability \ +and a higher incidence of recovery can give very close observed transition) it might \ +be very useful to look not only at linear confidence intervals estimated from the \ +variances but at the covariance matrix. And instead of looking at the estimated coefficients \ +(parameters) of the logistic regression, it might be more meaningful to visualize the \ +covariance matrix of the one-step probabilities. \ +See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres); - fprintf(fichtm," - Variance of one-step probabilities: %s
    \n", + fprintf(fichtm," - Standard deviation of one-step probabilities: %s
    \n", subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); fprintf(fichtm,"\ - Variance-covariance of one-step probabilities: %s
    \n", @@ -4356,26 +4619,26 @@ fprintf(fichtm," \n"); fflush(fichtm); @@ -4403,11 +4666,11 @@ void printinggnuplot(char fileres[], cha fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n"); for (cpt=1; cpt<= nlstate ; cpt ++) { for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ - 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,"\nset out \"%s%d_%d.svg\" \n",subdirf2(optionfilefiname,"v"),cpt,k1); + fprintf(ficgp,"\n#set out \"v%s%d_%d.svg\" \n",optionfilefiname,cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \n\ set ylabel \"Probability\" \n\ -set ter png small size 320, 240\n\ +set ter svg size 640, 480\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 ++) { @@ -4430,8 +4693,8 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u /*2 eme*/ fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); 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 size 320, 240\nplot [%.f:%.f] ",ageminpar,fage); + fprintf(ficgp,"\nset out \"%s%d.svg\" \n",subdirf2(optionfilefiname,"e"),k1); + fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); for (i=1; i<= nlstate+1 ; i ++) { k=2*i; @@ -4464,8 +4727,8 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u 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 size 320, 240\n\ + fprintf(ficgp,"\nset out \"%s%d%d.svg\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1); + fprintf(ficgp,"set ter svg size 640, 480\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) "); @@ -4489,9 +4752,9 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ k=3; fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt); - fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1); + fprintf(ficgp,"\nset out \"%s%d_%d.svg\" \n",subdirf2(optionfilefiname,"p"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ -set ter png small size 320, 240\n\ +set ter svg size 640, 480\n\ unset log y\n\ plot [%.f:%.f] ", ageminpar, agemaxpar); for (i=1; i<= nlstate ; i ++){ @@ -4510,56 +4773,93 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) } /* end covariate */ /* proba elementaires */ + fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); for(i=1,jk=1; i <=nlstate; i++){ + fprintf(ficgp,"# initial state %d\n",i); for(k=1; k <=(nlstate+ndeath); k++){ if (k != i) { + fprintf(ficgp,"# current state %d\n",k); for(j=1; j <=ncovmodel; j++){ - fprintf(ficgp,"p%d=%f ",jk,p[jk]); + fprintf(ficgp,"p%d=%f; ",jk,p[jk]); jk++; - fprintf(ficgp,"\n"); } + fprintf(ficgp,"\n"); } } } + fprintf(ficgp,"##############\n#\n"); + /*goto avoid;*/ + fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n"); + fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n"); + fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n"); + fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n"); + fprintf(ficgp,"# logi(p13/p11)=p6 +p7*age +p8*age*age+ p9*V1+ p10*V1*age\n"); + fprintf(ficgp,"# p12+p13+p14+p11=1=p11(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); + fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n"); + fprintf(ficgp,"# p11=1/(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); + fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n"); + fprintf(ficgp,"# p12=exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)/\n"); + fprintf(ficgp,"# (1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); + fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n"); + fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n"); + fprintf(ficgp,"#\n"); for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/ + fprintf(ficgp,"# ng=%d\n",ng); + fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m); for(jk=1; jk <=m; jk++) { - fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng); + fprintf(ficgp,"# jk=%d\n",jk); + fprintf(ficgp,"\nset out \"%s%d_%d.svg\" \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 size 320, 240\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar); + fprintf(ficgp,"\nset ter svg size 640, 480\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); + if(nagesqr==0) + fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); + else /* nagesqr =1 */ + fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr); else - fprintf(ficgp," exp(p%d+p%d*x",i,i+1); + if(nagesqr==0) + fprintf(ficgp," exp(p%d+p%d*x",i,i+1); + else /* nagesqr =1 */ + fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr); ij=1;/* To be checked else nbcode[0][0] wrong */ - for(j=3; j <=ncovmodel; j++) { - if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */ - fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); - ij++; + for(j=3; j <=ncovmodel-nagesqr; j++) { + /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ + if(ij <=cptcovage) { /* Bug valgrind */ + if((j-2)==Tage[ij]) { /* Bug valgrind */ + fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); + ij++; + } } else - fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]); + fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(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); + for(k1=1; k1 <=nlstate; k1++){ + if(nagesqr==0) + fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); + else /* nagesqr =1 */ + fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr); + 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++; + for(j=3; j <=ncovmodel-nagesqr; j++){ + if(ij <=cptcovage) { /* Bug valgrind */ + if((j-2)==Tage[ij]) { /* Bug valgrind */ + fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(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,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); } fprintf(ficgp,")"); } @@ -4686,7 +4986,7 @@ void prevforecast(char fileres[], double 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," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } fprintf(ficresf,"******\n"); fprintf(ficresf,"# Covariate valuofcovar yearproj age"); @@ -4710,7 +5010,7 @@ void prevforecast(char fileres[], double 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,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); } for(j=1; j<=nlstate+ndeath;j++) { @@ -4808,7 +5108,7 @@ void populforecast(char fileres[], doubl k=k+1; fprintf(ficrespop,"\n#******"); for(j=1;j<=cptcoveff;j++) { - fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } fprintf(ficrespop,"******\n"); fprintf(ficrespop,"# Age"); @@ -5123,7 +5423,7 @@ void printinghtmlmort(char fileres[], ch 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,"

    "); fprintf(fichtm,""); fprintf(fichtm,"