/* $Id$
$State$
$Log$
+ Revision 1.223 2016/02/19 09:23:35 brouard
+ Summary: temporary
+
Revision 1.222 2016/02/17 08:14:50 brouard
Summary: Probably last 0.98 stable version 0.98r6
/* #define DEBUGLINMIN */
/* #define DEBUGHESS */
#define DEBUGHESSIJ
-/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */
+/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan) *\/ */
#define POWELL /* Instead of NLOPT */
-#define POWELLF1F3 /* Skip test */
+#define POWELLNOF3INFF1TEST /* Skip test */
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
/* $State$ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-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 copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */
int cptcovprodnoage=0; /**< Number of covariate products without age */
int cptcoveff=0; /* Total number of covariates to vary for printing results */
+int ncoveff=0; /* Total number of effective covariates in the model */
+int nqveff=0; /**< nqveff number of effective quantitative variables */
+int ntveff=0; /**< ntveff number of effective time varying variables */
+int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */
int cptcov=0; /* Working variable */
int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
int npar=NPARMAX;
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]*age; */
+double **coqvar; /* Fixed quantitative covariate */
double ***cotvar; /* Time varying covariate */
double ***cotqvar; /* Time varying quantitative covariate */
-double **coqvar; /* Fixed quantitative covariate */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
int *Tage;
etemp=e;
e=d;
if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
- d=CGOLD*(e=(x >= xm ? a-x : b-x));
+ d=CGOLD*(e=(x >= xm ? a-x : b-x));
else {
- d=p/q;
- u=x+d;
- if (u-a < tol2 || b-u < tol2)
- d=SIGN(tol1,xm-x);
+ d=p/q;
+ u=x+d;
+ if (u-a < tol2 || b-u < tol2)
+ d=SIGN(tol1,xm-x);
}
} else {
d=CGOLD*(e=(x >= xm ? a-x : b-x));
} else {
if (u < x) a=u; else b=u;
if (fu <= fw || w == x) {
- v=w;
- w=u;
- fv=fw;
- fw=fu;
+ v=w;
+ w=u;
+ fv=fw;
+ fw=fu;
} else if (fu <= fv || v == x || v == w) {
- v=u;
- fv=fu;
+ v=u;
+ fv=fu;
}
}
}
*cx=(*bx)+GOLD*(*bx-*ax);
*fc=(*func)(*cx);
#ifdef DEBUG
- printf("mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);
- fprintf(ficlog,"mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);
+ printf("mnbrak0 a=%lf *fa=%lf, b=%lf *fb=%lf, c=%lf *fc=%lf\n",*ax,*fa,*bx,*fb,*cx, *fc);
+ fprintf(ficlog,"mnbrak0 a=%lf *fa=%lf, b=%lf *fb=%lf, c=%lf *fc=%lf\n",*ax,*fa,*bx,*fb,*cx, *fc);
#endif
- while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc */
+ while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc. If fc=inf it exits and if flat fb=fc it exits too.*/
r=(*bx-*ax)*(*fb-*fc);
- q=(*bx-*cx)*(*fb-*fa);
+ q=(*bx-*cx)*(*fb-*fa); /* What if fa=inf */
u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscissa of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */
ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscissa where function should be evaluated */
double A, fparabu;
A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u);
fparabu= *fa - A*(*ax-u)*(*ax-u);
- printf("mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
- fprintf(ficlog, "mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
+ printf("\nmnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f, q=%lf < %lf=r)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu,q,r);
+ fprintf(ficlog,"\nmnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f, q=%lf < %lf=r)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu,q,r);
/* And thus,it can be that fu > *fc even if fparabu < *fc */
/* mnbrak (*ax=7.666299858533, *fa=299039.693133272231), (*bx=8.595447774979, *fb=298976.598289369489),
(*cx=10.098840694817, *fc=298946.631474258087), (*u=9.852501168332, fu=298948.773013752128, fparabu=298945.434711494134) */
/* fu = *fc; */
/* *fc =dum; */
/* } */
-#ifdef DEBUG
- printf("mnbrak34 fu < or >= fc \n");
- fprintf(ficlog, "mnbrak34 fu < fc\n");
+#ifdef DEBUGMNBRAK
+ double A, fparabu;
+ A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u);
+ fparabu= *fa - A*(*ax-u)*(*ax-u);
+ printf("\nmnbrak35 ax=%lf fa=%lf bx=%lf fb=%lf, u=%lf fp=%lf fu=%lf < or >= fc=%lf cx=%lf, q=%lf < %lf=r \n",*ax, *fa, *bx,*fb,u,fparabu,fu,*fc,*cx,q,r);
+ fprintf(ficlog,"\nmnbrak35 ax=%lf fa=%lf bx=%lf fb=%lf, u=%lf fp=%lf fu=%lf < or >= fc=%lf cx=%lf, q=%lf < %lf=r \n",*ax, *fa, *bx,*fb,u,fparabu,fu,*fc,*cx,q,r);
#endif
dum=u; /* Shifting c and u */
u = *cx;
#endif
} else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */
#ifdef DEBUG
- printf("mnbrak2 u after c but before ulim\n");
- fprintf(ficlog, "mnbrak2 u after c but before ulim\n");
+ printf("\nmnbrak2 u=%lf after c=%lf but before ulim\n",u,*cx);
+ fprintf(ficlog,"\nmnbrak2 u=%lf after c=%lf but before ulim\n",u,*cx);
#endif
fu=(*func)(u);
if (fu < *fc) {
#ifdef DEBUG
- printf("mnbrak2 u after c but before ulim AND fu < fc\n");
- fprintf(ficlog, "mnbrak2 u after c but before ulim AND fu <fc \n");
+ printf("\nmnbrak2 u=%lf after c=%lf but before ulim=%lf AND fu=%lf < %lf=fc\n",u,*cx,ulim,fu, *fc);
+ fprintf(ficlog,"\nmnbrak2 u=%lf after c=%lf but before ulim=%lf AND fu=%lf < %lf=fc\n",u,*cx,ulim,fu, *fc);
+#endif
+ SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
+ SHFT(*fb,*fc,fu,(*func)(u))
+#ifdef DEBUG
+ printf("\nmnbrak2 shift GOLD c=%lf",*cx+GOLD*(*cx-*bx));
#endif
- SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
- SHFT(*fb,*fc,fu,(*func)(u))
}
} else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */
#ifdef DEBUG
- printf("mnbrak2 u outside ulim (verifying that ulim is beyond c)\n");
- fprintf(ficlog, "mnbrak2 u outside ulim (verifying that ulim is beyond c)\n");
+ printf("\nmnbrak2 u=%lf outside ulim=%lf (verifying that ulim is beyond c=%lf)\n",u,ulim,*cx);
+ fprintf(ficlog,"\nmnbrak2 u=%lf outside ulim=%lf (verifying that ulim is beyond c=%lf)\n",u,ulim,*cx);
#endif
u=ulim;
fu=(*func)(u);
} else { /* u could be left to b (if r > q parabola has a maximum) */
#ifdef DEBUG
- printf("mnbrak2 u could be left to b (if r > q parabola has a maximum)\n");
- fprintf(ficlog, "mnbrak2 u could be left to b (if r > q parabola has a maximum)\n");
+ printf("\nmnbrak2 u=%lf could be left to b=%lf (if r=%lf > q=%lf parabola has a maximum)\n",u,*bx,r,q);
+ fprintf(ficlog,"\nmnbrak2 u=%lf could be left to b=%lf (if r=%lf > q=%lf parabola has a maximum)\n",u,*bx,r,q);
#endif
u=(*cx)+GOLD*(*cx-*bx);
fu=(*func)(u);
+#ifdef DEBUG
+ printf("\nmnbrak2 new u=%lf fu=%lf shifted gold left from c=%lf and b=%lf \n",u,fu,*cx,*bx);
+ fprintf(ficlog,"\nmnbrak2 new u=%lf fu=%lf shifted gold left from c=%lf and b=%lf \n",u,fu,*cx,*bx);
+#endif
} /* end tests */
SHFT(*ax,*bx,*cx,u)
SHFT(*fa,*fb,*fc,fu)
#ifdef DEBUG
- printf("mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu);
- fprintf(ficlog, "mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu);
+ printf("\nmnbrak2 shift (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc);
+ fprintf(ficlog, "\nmnbrak2 shift (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc);
#endif
} /* end while; ie return (a, b, c, fa, fb, fc) such that a < b < c with f(a) > f(b) and fb < f(c) */
}
double *pcom,*xicom;
double (*nrfunc)(double []);
+#ifdef LINMINORIGINAL
void linmin(double p[], double xi[], int n, double *fret,double (*func)(double []))
+#else
+void linmin(double p[], double xi[], int n, double *fret,double (*func)(double []), int *flat)
+#endif
{
double brent(double ax, double bx, double cx,
double (*f)(double), double tol, double *xmin);
#ifdef LINMINORIGINAL
#else
if (fx != fx){
- xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */
- printf("|");
- fprintf(ficlog,"|");
+ 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);
+ 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);
+ }while(fx != fx && xxs > 1.e-5);
#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);
#endif
+#ifdef LINMINORIGINAL
+#else
+ if(fb == fx){ /* Flat function in the direction */
+ xmin=xx;
+ *flat=1;
+ }else{
+ *flat=0;
+#endif
+ /*Flat mnbrak2 shift (*ax=0.000000000000, *fa=51626.272983130431), (*bx=-1.618034000000, *fb=51590.149499362531), (*cx=-4.236068025156, *fc=51590.149499362531) */
*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);
+ printf("retour brent from bracket (a=%lf fa=%lf, xx=%lf fx=%lf, b=%lf fb=%lf): fret=%lf xmin=%lf\n",ax,fa,xx,fx,bx,fb,*fret,xmin);
+ fprintf(ficlog,"retour brent from bracket (a=%lf fa=%lf, xx=%lf fx=%lf, b=%lf fb=%lf): fret=%lf xmin=%lf\n",ax,fa,xx,fx,bx,fb,*fret,xmin);
+#endif
+#ifdef LINMINORIGINAL
+#else
+ }
#endif
#ifdef DEBUGLINMIN
printf("linmin end ");
output, p is set to the best point found, xi is the then-current direction set, fret is the returned
function value at p , and iter is the number of iterations taken. The routine linmin is used.
*/
+#ifdef LINMINORIGINAL
+#else
+ int *flatdir; /* Function is vanishing in that direction */
+ int flat=0; /* Function is vanishing in that direction */
+#endif
void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,
double (*func)(double []))
{
- void linmin(double p[], double xi[], int n, double *fret,
+#ifdef LINMINORIGINAL
+ void linmin(double p[], double xi[], int n, double *fret,
double (*func)(double []));
+#else
+ void linmin(double p[], double xi[], int n, double *fret,
+ double (*func)(double []),int *flat);
+#endif
int i,ibig,j;
double del,t,*pt,*ptt,*xit;
double directest;
double fp,fptt;
double *xits;
int niterf, itmp;
+#ifdef LINMINORIGINAL
+#else
+
+ flatdir=ivector(1,n);
+ for (j=1;j<=n;j++) flatdir[j]=0;
+#endif
pt=vector(1,n);
ptt=vector(1,n);
rforecast_time=rcurr_time;
itmp = strlen(strcurr);
if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */
- strcurr[itmp-1]='\0';
+ strcurr[itmp-1]='\0';
printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
for(niterf=10;niterf<=30;niterf+=10){
- rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
- forecast_time = *localtime(&rforecast_time);
- strcpy(strfor,asctime(&forecast_time));
- itmp = strlen(strfor);
- if(strfor[itmp-1]=='\n')
- strfor[itmp-1]='\0';
- printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
- fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+ rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
+ forecast_time = *localtime(&rforecast_time);
+ strcpy(strfor,asctime(&forecast_time));
+ itmp = strlen(strfor);
+ if(strfor[itmp-1]=='\n')
+ strfor[itmp-1]='\0';
+ printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+ fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
}
}
for (i=1;i<=n;i++) { /* For each direction i */
#endif
printf("%d",i);fflush(stdout); /* print direction (parameter) i */
fprintf(ficlog,"%d",i);fflush(ficlog);
+#ifdef LINMINORIGINAL
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 */
+#else
+ linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
+ flatdir[i]=flat; /* Function is vanishing in that direction i */
+#endif
+ /* 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;
+ /* 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;
}
#ifdef DEBUG
printf("%d %.12e",i,(*fret));
fprintf(ficlog,"%d %.12e",i,(*fret));
for (j=1;j<=n;j++) {
- xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
- printf(" x(%d)=%.12e",j,xit[j]);
- fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
+ xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
+ printf(" x(%d)=%.12e",j,xit[j]);
+ fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
}
for(j=1;j<=n;j++) {
- printf(" p(%d)=%.12e",j,p[j]);
- fprintf(ficlog," p(%d)=%.12e",j,p[j]);
+ printf(" p(%d)=%lf ",j,p[j]);
+ fprintf(ficlog," p(%d)=%lf ",j,p[j]);
}
printf("\n");
fprintf(ficlog,"\n");
/* 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) */
+ for(j=1;j<=n;j++) {
+ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
+
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 */
/* 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 */
}
#endif
-
+#ifdef LINMINORIGINAL
+#else
+ free_ivector(flatdir,1,n);
+#endif
free_vector(xit,1,n);
free_vector(xits,1,n);
free_vector(ptt,1,n);
pt[j]=p[j];
}
fptt=(*func)(ptt); /* f_3 */
-#ifdef POWELLF1F3
+#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */
+ if (*iter <=4) {
+#else
+#ifdef POWELLNOF3INFF1TEST /* skips test F3 <F1 */
#else
if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
#endif
/* Let f"(x2) be the 2nd derivative equal everywhere. */
/* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */
/* will reach at f3 = fm + h^2/2 f"m ; f" = (f1 -2f2 +f3 ) / h**2 */
- /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */
+ /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del or directest <0 */
+ /* also lamda^2=(f1-f2)^2/mu² is a parasite solution of powell */
+ /* For powell, inclusion of this average direction is only if t(del)<0 or del inbetween mu^2 and lambda^2 */
/* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
+ /* Even if f3 <f1, directest can be negative and t >0 */
+ /* mu² and del² are equal when f3=f1 */
+ /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
+ /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
+ /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */
+ /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */
#ifdef NRCORIGINAL
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
#else
if (t < 0.0) { /* Then we use it for new direction */
#else
if (directest*t < 0.0) { /* Contradiction between both tests */
- printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
+ printf("directest= %.12lf (if <0 we include P0 Pn as new direction), 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 (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
+ fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), 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
#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]);
- fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0){
- printf("\n");
- fprintf(ficlog,"\n");
- }
- }
+ 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]);
+ fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
+ printf("\n");
+ fprintf(ficlog,"\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]);
- fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0){
- printf("\n");
- fprintf(ficlog,"\n");
- }
- }
+#ifdef LINMINORIGINAL
+ linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+#else
+ linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+ flatdir[i]=flat; /* Function is vanishing in that direction i */
#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 */
- }
- printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
- fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+#ifdef DEBUGLINMIN
+ for (j=1;j<=n;j++) {
+ printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
+ printf("\n");
+ fprintf(ficlog,"\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 */
+ }
+#ifdef LINMINORIGINAL
+#else
+ printf("Flat directions\n");
+ fprintf(ficlog,"Flat directions\n");
+ for (j=1;j<=n;j++) {
+ printf("flatdir[%d]=%d ",j,flatdir[j]);
+ fprintf(ficlog,"flatdir[%d]=%d ",j,flatdir[j]);
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
+#endif
+ printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+ fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+
#ifdef DEBUG
- printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
- fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
- for(j=1;j<=n;j++){
- printf(" %.12e",xit[j]);
- fprintf(ficlog," %.12e",xit[j]);
- }
- printf("\n");
- fprintf(ficlog,"\n");
+ printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+ fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+ for(j=1;j<=n;j++){
+ printf(" %lf",xit[j]);
+ fprintf(ficlog," %lf",xit[j]);
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
#endif
} /* end of t or directest negative */
-#ifdef POWELLF1F3
+#ifdef POWELLNOF3INFF1TEST
#else
} /* end if (fptt < fp) */
+#endif
+ } /*NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done */
#endif
} /* loop iteration */
}
/*************** log-likelihood *************/
double func( double *x)
{
- int i, ii, j, k, mi, d, kk;
- int ioffset;
- double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
- double **out;
- double sw; /* Sum of weights */
- double lli; /* Individual log likelihood */
- int s1, s2;
- int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate */
- 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]);*/
- /*for(i=1;i<imx;i++)
- printf(" %d\n",s[4][i]);
- */
-
- ++countcallfunc;
-
- cov[1]=1.;
-
- for(k=1; k<=nlstate; k++) ll[k]=0.;
+ int i, ii, j, k, mi, d, kk;
+ int ioffset=0;
+ double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
+ double **out;
+ double sw; /* Sum of weights */
+ double lli; /* Individual log likelihood */
+ int s1, s2;
+ int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate */
+ 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]);*/
+ /*for(i=1;i<imx;i++)
+ printf(" %d\n",s[4][i]);
+ */
- if(mle==1){
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- /* Computes the values of the ncovmodel covariates of the model
- depending if the covariates are fixed or varying (age dependent) and stores them in cov[]
- Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
- to be observed in j being in i according to the model.
- */
- ioffset=2+nagesqr;
- for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */
- cov[++ioffset]=covar[Tvar[k]][i];
- }
- for(iqv=1; iqv <= nqv; iqv++){ /* Varying quantitatives covariates */
- /* cov[2+nagesqr+cptcovn+iqv]=varq[mw[mi+1][i]][iqv][i]; */
- }
+ ++countcallfunc;
+
+ cov[1]=1.;
+
+ for(k=1; k<=nlstate; k++) ll[k]=0.;
+ ioffset=0;
+ if(mle==1){
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ /* Computes the values of the ncovmodel covariates of the model
+ depending if the covariates are fixed or varying (age dependent) and stores them in cov[]
+ Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
+ to be observed in j being in i according to the model.
+ */
+ ioffset=2+nagesqr+cptcovage;
+ /* for (k=1; k<=cptcovn;k++){ /\* Simple and product covariates without age* products *\/ */
+ for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */
+ cov[++ioffset]=covar[Tvar[k]][i];
+ }
+ for(iqv=1; iqv <= nqveff; iqv++){ /* Quantitatives covariates */
+ cov[++ioffset]=coqvar[iqv][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]
- has been calculated etc */
- /* For an individual i, wav[i] gives the number of effective waves */
- /* We compute the contribution to Likelihood of each effective transition
- mw[mi][i] is real wave of the mi th effectve wave */
- /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
- s2=s[mw[mi+1][i]][i];
- And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]
- But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
- meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
- */
- for(mi=1; mi<= wav[i]-1; mi++){
- for(itv=1; itv <= ntv; itv++){ /* Varying dummy covariates */
- cov[++ioffset]=cotvar[mw[mi+1][i]][itv][i];
- }
- for(iqtv=1; iqtv <= nqtv; iqtv++){ /* Varying quantitatives covariates */
- /* cov[2+nagesqr+cptcovn+nqv+ntv+iqtv]=varq[mw[mi+1][i]][iqtv][i]; */
- }
- ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv;
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact; /* Should be changed here */
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
- }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
-
- /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
- /* But now since version 0.9 we anticipate for bias at large stepm.
- * If stepm is larger than one month (smallest stepm) and if the exact delay
- * (in months) between two waves is not a multiple of stepm, we rounded to
- * the nearest (and in case of equal distance, to the lowest) interval but now
- * we keep into memory the bias bh[mi][i] and also the previous matrix product
- * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
- * probability in order to take into account the bias as a fraction of the way
- * from savm to out if bh is negative or even beyond if bh is positive. bh varies
- * -stepm/2 to stepm/2 .
- * For stepm=1 the results are the same as for previous versions of Imach.
- * For stepm > 1 the results are less biased than in previous versions.
- */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- bbh=(double)bh[mi][i]/(double)stepm;
- /* bias bh is positive if real duration
- * is higher than the multiple of stepm and negative otherwise.
- */
- /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
- if( s2 > nlstate){
- /* i.e. if s2 is a death state and if the date of death is known
- then the contribution to the likelihood is the probability to
- die between last step unit time and current step unit time,
- which is also equal to probability to die before dh
- minus probability to die before dh-stepm .
- In version up to 0.92 likelihood was computed
- as if date of death was unknown. Death was treated as any other
- health state: the date of the interview describes the actual state
- and not the date of a change in health state. The former idea was
- to consider that at each interview the state was recorded
- (healthy, disable or death) and IMaCh was corrected; but when we
- introduced the exact date of death then we should have modified
- the contribution of an exact death to the likelihood. This new
- contribution is smaller and very dependent of the step unit
- stepm. It is no more the probability to die between last interview
- and month of death but the probability to survive from last
- interview up to one month before death multiplied by the
- probability to die within a month. Thanks to Chris
- Jackson for correcting this bug. Former versions increased
- mortality artificially. The bad side is that we add another loop
- which slows down the processing. The difference can be up to 10%
- lower mortality.
- */
- /* If, at the beginning of the maximization mostly, the
- cumulative probability or probability to be dead is
- constant (ie = 1) over time d, the difference is equal to
- 0. out[s1][3] = savm[s1][3]: probability, being at state
- s1 at precedent wave, to be dead a month before current
- wave is equal to probability, being at state s1 at
- precedent wave, to be dead at mont of the current
- wave. Then the observed probability (that this person died)
- is null according to current estimated parameter. In fact,
- it should be very low but not zero otherwise the log go to
- infinity.
- */
+ /* 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]
+ has been calculated etc */
+ /* For an individual i, wav[i] gives the number of effective waves */
+ /* We compute the contribution to Likelihood of each effective transition
+ mw[mi][i] is real wave of the mi th effectve wave */
+ /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
+ s2=s[mw[mi+1][i]][i];
+ And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]
+ But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
+ meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
+ */
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */
+ cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i];
+ }
+ for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */
+ cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i];
+ }
+ /* ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; */
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact; /* Should be changed here */
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
+
+ /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+ /* But now since version 0.9 we anticipate for bias at large stepm.
+ * If stepm is larger than one month (smallest stepm) and if the exact delay
+ * (in months) between two waves is not a multiple of stepm, we rounded to
+ * the nearest (and in case of equal distance, to the lowest) interval but now
+ * we keep into memory the bias bh[mi][i] and also the previous matrix product
+ * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
+ * probability in order to take into account the bias as a fraction of the way
+ * from savm to out if bh is negative or even beyond if bh is positive. bh varies
+ * -stepm/2 to stepm/2 .
+ * For stepm=1 the results are the same as for previous versions of Imach.
+ * For stepm > 1 the results are less biased than in previous versions.
+ */
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ /* bias bh is positive if real duration
+ * is higher than the multiple of stepm and negative otherwise.
+ */
+ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
+ if( s2 > nlstate){
+ /* i.e. if s2 is a death state and if the date of death is known
+ then the contribution to the likelihood is the probability to
+ die between last step unit time and current step unit time,
+ which is also equal to probability to die before dh
+ minus probability to die before dh-stepm .
+ In version up to 0.92 likelihood was computed
+ as if date of death was unknown. Death was treated as any other
+ health state: the date of the interview describes the actual state
+ and not the date of a change in health state. The former idea was
+ to consider that at each interview the state was recorded
+ (healthy, disable or death) and IMaCh was corrected; but when we
+ introduced the exact date of death then we should have modified
+ the contribution of an exact death to the likelihood. This new
+ contribution is smaller and very dependent of the step unit
+ stepm. It is no more the probability to die between last interview
+ and month of death but the probability to survive from last
+ interview up to one month before death multiplied by the
+ probability to die within a month. Thanks to Chris
+ Jackson for correcting this bug. Former versions increased
+ mortality artificially. The bad side is that we add another loop
+ which slows down the processing. The difference can be up to 10%
+ lower mortality.
+ */
+ /* If, at the beginning of the maximization mostly, the
+ cumulative probability or probability to be dead is
+ constant (ie = 1) over time d, the difference is equal to
+ 0. out[s1][3] = savm[s1][3]: probability, being at state
+ s1 at precedent wave, to be dead a month before current
+ wave is equal to probability, being at state s1 at
+ precedent wave, to be dead at mont of the current
+ wave. Then the observed probability (that this person died)
+ is null according to current estimated parameter. In fact,
+ it should be very low but not zero otherwise the log go to
+ infinity.
+ */
/* #ifdef INFINITYORIGINAL */
/* lli=log(out[s1][s2] - savm[s1][s2]); */
/* #else */
/* else */
/* lli=log(out[s1][s2] - savm[s1][s2]); */
/* #endif */
- lli=log(out[s1][s2] - savm[s1][s2]);
+ 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 */
- }
- /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
- /*if(lli ==000.0)*/
- /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- /* if (lli < log(mytinydouble)){ */
- /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */
- /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */
- /* } */
- } /* end of wave */
- } /* end of individual */
- } else if(mle==2){
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<=dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ } 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 */
+ }
+ /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
+ /*if(lli ==000.0)*/
+ /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ /* if (lli < log(mytinydouble)){ */
+ /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */
+ /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */
+ /* } */
+ } /* end of wave */
+ } /* end of individual */
+ } else if(mle==2){
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<=dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- bbh=(double)bh[mi][i]/(double)stepm;
- lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- } /* end of wave */
- } /* end of individual */
- } else if(mle==3){ /* exponential inter-extrapolation */
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ } /* end of wave */
+ } /* end of individual */
+ } else if(mle==3){ /* exponential inter-extrapolation */
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- bbh=(double)bh[mi][i]/(double)stepm;
- lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- } /* end of wave */
- } /* end of individual */
- }else if (mle==4){ /* ml=4 no inter-extrapolation */
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ } /* end of wave */
+ } /* end of individual */
+ }else if (mle==4){ /* ml=4 no inter-extrapolation */
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- 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 */
- }
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ s1=s[mw[mi][i]][i];
+ 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 */
+ }
+ ipmx +=1;
+ sw += weight[i];
+ 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]); */
- } /* end of wave */
- } /* end of individual */
- }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
+ } /* end of wave */
+ } /* end of individual */
+ }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
- ipmx +=1;
- sw += weight[i];
- 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]);*/
- } /* end of wave */
- } /* end of individual */
- } /* End of if */
- for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
- /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
- l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
- return -l;
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
+ ipmx +=1;
+ sw += weight[i];
+ 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]);*/
+ } /* end of wave */
+ } /* end of individual */
+ } /* End of if */
+ for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
+ /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
+ l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+ return -l;
}
/*************** log-likelihood *************/
{
/* Same as likeli but slower because of a lot of printf and if */
int i, ii, j, k, mi, d, kk;
+ int ioffset=0;
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
double lli; /* Individual log likelihood */
double llt;
int s1, s2;
+ int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate */
double bbh, survp;
double agexact;
double agebegin, ageend;
cov[1]=1.;
for(k=1; k<=nlstate; k++) ll[k]=0.;
-
+ ioffset=0;
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ ioffset=2+nagesqr+cptcovage;
+ /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
+ for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */
+ cov[++ioffset]=covar[Tvar[k]][i];
+ }
+ for(iqv=1; iqv <= nqveff; iqv++){ /* Quantitatives covariates */
+ cov[++ioffset]=coqvar[iqv][i];
+ }
+
for(mi=1; mi<= wav[i]-1; mi++){
+ for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */
+ cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i];
+ }
+ for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */
+ cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i];
+ }
for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
for(d=0; d<dh[mi][i]; d++){ /* Delay between two effective waves */
- /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
- and mw[mi+1][i]. dh depends on stepm.*/
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
-
- /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */
- /* 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */
- savm=oldm;
- oldm=newm;
+ /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
+ and mw[mi+1][i]. dh depends on stepm.*/
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
+
+ /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */
+ /* 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */
+ savm=oldm;
+ oldm=newm;
} /* end mult */
s1=s[mw[mi][i]][i];
* is higher than the multiple of stepm and negative otherwise.
*/
if( s2 > nlstate && (mle <5) ){ /* Jackson */
- lli=log(out[s1][s2] - savm[s1][s2]);
+ 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];
- lli= log(survp);
+ for (j=1,survp=0. ; j<=nlstate; j++)
+ survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ lli= log(survp);
}else if (mle==1){
- lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
} else if(mle==2){
- lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
+ lli= (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 */
} else if(mle==3){ /* exponential inter-extrapolation */
- lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
+ lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
} else if (mle==4){ /* mle=4 no inter-extrapolation */
- lli=log(out[s1][s2]); /* Original formula */
+ lli=log(out[s1][s2]); /* Original formula */
} else{ /* mle=0 back to 1 */
- lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
- /*lli=log(out[s1][s2]); */ /* Original formula */
+ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+ /*lli=log(out[s1][s2]); */ /* Original formula */
} /* End of if */
ipmx +=1;
sw += weight[i];
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 %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
+ fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
%11.6f %11.6f %11.6f ", \
- 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;
- fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
- }
- fprintf(ficresilk," %10.6f\n", -llt);
+ 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;
+ fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
+ }
+ fprintf(ficresilk," %10.6f\n", -llt);
}
} /* end of wave */
} /* end of individual */
int mi; /* Effective wave */
int first;
double ***freq; /* Frequencies */
+ double *meanq;
+ double **meanqt;
double *pp, **prop, *posprop, *pospropt;
double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0;
char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];
posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */
pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */
/* prop=matrix(1,nlstate,iagemin,iagemax+3); */
+ meanq=vector(1,nqveff);
+ meanqt=matrix(1,lastpass,1,nqtveff);
strcpy(fileresp,"P_");
strcat(fileresp,fileresu);
/*strcat(fileresphtm,fileresu);*/
freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
j1=0;
- j=cptcoveff;
+ j=ncoveff;
if (cptcovn<1) {j=1;ncodemax[1]=1;}
first=1;
Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff
*/
- for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */
+ for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination excluding varying and quantitatives */
posproptt=0.;
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/
posprop[i]=0;
pospropt[i]=0;
}
+ for (z1=1; z1<= nqveff; z1++) {
+ meanq[z1]+=0.;
+ for(m=1;m<=lastpass;m++){
+ meanqt[m][z1]=0.;
+ }
+ }
dateintsum=0;
k2cpt=0;
-
+ /* For that comination of covariate j1, we count and print the frequencies */
for (iind=1; iind<=imx; iind++) { /* For each individual iind */
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 (nqveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+ for (z1=1; z1<= nqveff; z1++) {
+ meanq[z1]+=coqvar[Tvar[z1]][iind];
+ }
+ for (z1=1; z1<=ncoveff; z1++) {
+ /* if(Tvaraff[z1] ==-20){ */
+ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
+ /* }else if(Tvaraff[z1] ==-10){ */
+ /* /\* sumnew+=coqvar[z1][iind]; *\/ */
+ /* }else */
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
- /* Tests if the value of each of the covariates of i is equal to filter j1 */
+ /* Tests if this individual i responded to j1 (V4=1 V3=0) */
bool=0;
/* 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),
} /* end z1 */
} /* cptcovn > 0 */
- if (bool==1){
+ if (bool==1){ /* We selected an individual iin satisfying combination j1 */
/* for(m=firstpass; m<=lastpass; m++){ */
for(mi=1; mi<wav[iind];mi++){
m=mw[mi][iind];
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
pstamp(ficresp);
- if (cptcovn>0) {
+ if (ncoveff>0) {
fprintf(ficresp, "\n#********** Variable ");
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++){
+ for (z1=1; z1<=ncoveff; 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(ficresphtm, "**********</h3>\n");
fprintf(ficresphtmfr, "**********</h3>\n");
fprintf(ficlog, "\n#********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ for (z1=1; z1<=ncoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficlog, "**********\n");
}
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
fclose(ficresp);
fclose(ficresphtm);
fclose(ficresphtmfr);
+ free_vector(meanq,1,nqveff);
+ free_matrix(meanqt,1,lastpass,1,nqtveff);
free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE);
free_vector(pospropt,1,nlstate);
free_vector(posprop,1,nlstate);
free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);
free_vector(pp,1,nlstate);
- /* End of Freq */
+ /* End of freqsummary */
}
/************ Prevalence ********************/
if (cptcovn<1) {j=1;ncodemax[1]=1;}
first=1;
- for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
+ for(j1=1; j1<= (int) pow(2,nqveff);j1++){ /* For each combination of covariate */
for (i=1; i<=nlstate; i++)
for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; 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++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
+ for (z1=1; z1<=nqveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
bool=0;
}
and mw[mi+1][i]. dh depends on stepm.
*/
- int i, mi, m;
+ int i=0, mi=0, m=0, mli=0;
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
double sum=0., jmean=0.;*/
- int first, firstwo, firsthree, firstfour;
+ int first=0, firstwo=0, firsthree=0, firstfour=0, firstfiv=0;
int j, k=0,jk, ju, jl;
double sum=0.;
first=0;
jmin=100000;
jmax=-1;
jmean=0.;
+
+/* Treating live states */
for(i=1; i<=imx; i++){ /* For simple cases and if state is death */
- mi=0;
+ mi=0; /* First valid wave */
+ mli=0; /* Last valid wave */
m=firstpass;
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 */
+ if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */
+ mli=m-1;/* mw[++mi][i]=m-1; */
+ }else 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;
+ mli=m;
+ } /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */
+ if(m < lastpass){ /* m < lastpass, standard case */
+ m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */
}
- if(m >=lastpass){
+ else{ /* m >= lastpass, eventual special issue with warning */
+#ifdef UNKNOWNSTATUSNOTCONTRIBUTING
+ break;
+#else
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);
+ printf("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 as pi. .\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;
}
- 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.\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 as pi. .\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;
+ mli=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);
+ 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.\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++;
+#endif
+ }/* End m >= lastpass */
}/* end while */
-
+
+ /* mi is the last effective wave, m is lastpass, mw[j][i] gives the # of j-th effective wave for individual i */
/* After last pass */
+/* Treating death states */
if (s[m][i] > nlstate){ /* In a death state */
+ /* if( mint[m][i]==mdc[m][i] && anint[m][i]==andc[m][i]){ /\* same date of death and date of interview *\/ */
+ /* } */
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 */
+ }
+#ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE
+ else if ((int) andc[i] != 9999) { /* Status is 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 ((int)anint[m][i]!= 9999) { /* date of last interview is known */
- 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 );
- firstwo=1;
+ if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */
+ nbwarn++;
+ if(firstfiv==0){
+ printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing 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], s[m][i], i,m );
+ firstfiv=1;
+ }else{
+ fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m );
+ }
+ }else{ /* Death occured afer last wave potential bias */
+ nberr++;
+ if(firstwo==0){
+ printf("Error! Death for individual %ld line=%d occurred at %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;
+ }
+ fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %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.\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.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
}else{ /* end date of interview is known */
/* death is known but not confirmed by death status at any wave */
if(firstfour==0){
}
fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %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.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
}
- }
- wav[i]=mi;
+ } /* end if date of death is known */
+#endif
+ wav[i]=mi; /* mi should be the last effective wave (or mli) */
+ /* wav[i]=mw[mi][i]; */
if(mi==0){
nbwarn++;
if(first==0){
/**< Uses cptcovn+2*cptcovprod as the number of covariates */
/* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1
* Boring subroutine which should only output nbcode[Tvar[j]][k]
- * Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2)
- * nbcode[Tvar[5]][1]= nbcode[2][1]=0, nbcode[2][2]=1 (usually);
+ * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable
+ * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually);
*/
int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
/* cptcoveff=0; */
- *cptcov=0;
+ /* *cptcov=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 */
+ /* Loop on covariates without age and products and no quantitative variable */
+ /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
+ for (j=1; j<=(*cptcov); 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
- * If product of Vn*Vm, still boolean *:
- * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
- * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */
- /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
- modality of the nth covariate of individual i. */
- if (ij > modmaxcovj)
- modmaxcovj=ij;
- else if (ij < modmincovj)
- modmincovj=ij;
- if ((ij < -1) && (ij > NCOVMAX)){
- printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
- exit(1);
- }else
- Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
+ modality of this covariate Vj*/
+ if(Tvar[j] >=1 && Tvar[j] <= *cptcov){ /* A real fixed covariate */
+ ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
+ * If product of Vn*Vm, still boolean *:
+ * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
+ * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */
+ /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
+ modality of the nth covariate of individual i. */
+ if (ij > modmaxcovj)
+ modmaxcovj=ij;
+ else if (ij < modmincovj)
+ modmincovj=ij;
+ if ((ij < -1) && (ij > NCOVMAX)){
+ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
+ exit(1);
+ }else
+ Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
/* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
/*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
/* 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 */
+ (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
+ female ies 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=0; i<=cptcode; i++) {*/
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]);
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 */
+ * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
} /* Ndum[-1] number of undefined modalities */
/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
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 */
+ 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 */
/* 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]++; /* Might be supersed V1 + V1*age */
- }
+ } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 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) */
- }else{
- /* Tvaraff[ij]=0; */
+ Tvaraff[++ij]=i; /*For printing (unclear) */
+ }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){
+ Tvaraff[++ij]=-10; /* Dont'n know how to treat quantitative variables yet */
+ }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){
+ Tvaraff[++ij]=i; /*For printing (unclear) */
+ }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){
+ Tvaraff[++ij]=-20; /* Dont'n know how to treat quantitative variables yet */
}
- }
+ } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
/* ij--; */
/* cptcoveff=ij; /\*Number of total covariates*\/ */
- *cptcov=ij; /*Number of total covariates*/
-
+ *cptcov=ij; /*Number of total real effective covariates: effective
+ * because they can be excluded from the model and real
+ * if in the model but excluded because missing values*/
}
cov[1]=1;
/* tj=cptcoveff; */
- tj = (int) pow(2,cptcoveff);
+ tj = (int) pow(2,nqveff);
if (cptcovn<1) {tj=1;ncodemax[1]=1;}
j1=0;
- for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates */
+ for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates or only once*/
if (cptcovn>0) {
fprintf(ficresprob, "\n#********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ for (z1=1; z1<=nqveff; 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]][codtabm(j1,z1)]);
+ for (z1=1; z1<=nqveff; 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]][codtabm(j1,z1)]);
+ for (z1=1; z1<=nqveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficgp, "**********\n#\n");
fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ for (z1=1; z1<=nqveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
fprintf(ficresprobcor, "\n#********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ for (z1=1; z1<=nqveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficresprobcor, "**********\n#");
if(invalidvarcomb[j1]){
fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1);
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
- m=pow(2,cptcoveff);
+ m=pow(2,nqveff);
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
jj1=0;
jj1++;
if (cptcovn > 0) {
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
- for (cpt=1; cpt<=cptcoveff;cpt++){
+ for (cpt=1; cpt<=nqveff;cpt++){
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout);
}
fflush(fichtm);
fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
- m=pow(2,cptcoveff);
+ m=pow(2,nqveff);
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
jj1=0;
jj1++;
if (cptcovn > 0) {
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
- for (cpt=1; cpt<=cptcoveff;cpt++)
+ for (cpt=1; cpt<=nqveff;cpt++)
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
/*#ifdef windows */
fprintf(ficgp,"cd \"%s\" \n",pathc);
/*#endif */
- m=pow(2,cptcoveff);
+ m=pow(2,nqveff);
/* Contribution to likelihood */
/* Plot the probability implied in the likelihood */
for (k1=1; k1<= m ; k1 ++) { /* For each valid combination of covariate */
/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
- for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
+ for (k=1; k<=nqveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */
+ lv= decodtabm(k1,k,nqveff); /* Should be the value of the covariate corresponding to k1 combination */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
- if(cptcoveff ==0){
+ if(nqveff ==0){
fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ", 2+(cpt-1), cpt );
}else{
kl=0;
- for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+ for (k=1; k<=nqveff; k++){ /* For each combination of covariate */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
/*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
- if(k==cptcoveff){
+ if(k==nqveff){
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' with line ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
6+(cpt-1), cpt );
}else{
for (k1=1; k1<= m ; k1 ++) {
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ for (k=1; k<=nqveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate number corresponding to k1 combination */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
for (cpt=1; cpt<= nlstate ; cpt ++) {
fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ for (k=1; k<=nqveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate number corresponding to k1 combination */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ for (k=1; k<=nqveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate number corresponding to k1 combination */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */
fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ for (k=1; k<=nqveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate number corresponding to k1 combination */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ for (k=1; k<=nqveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate number corresponding to k1 combination */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ for (k=1; k<=nqveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate number corresponding to k1 combination */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+ for (k=1; k<=nqveff; k++){ /* For each correspondig covariate value */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
}else{
fprintf(ficgp,",\\\n '' ");
}
- if(cptcoveff ==0){ /* No covariate */
+ if(nqveff ==0){ /* No covariate */
ioffset=2; /* Age is in 2 */
/*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
}else{ /* more than 2 covariates */
- if(cptcoveff ==1){
+ if(nqveff ==1){
ioffset=4; /* Age is in 4 */
}else{
ioffset=6; /* Age is in 6 */
fprintf(ficgp," u %d:(",ioffset);
kl=0;
strcpy(gplotcondition,"(");
- for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+ for (k=1; k<=nqveff; k++){ /* For each covariate writing the chain of conditions */
+ lv= decodtabm(k1,k,nqveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
kl++;
sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
kl++;
- if(k <cptcoveff && cptcoveff>1)
+ if(k <nqveff && nqveff>1)
sprintf(gplotcondition+strlen(gplotcondition)," && ");
}
strcpy(gplotcondition+strlen(gplotcondition),")");
fprintf(ficgp,"#\n");
for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
fprintf(ficgp,"# ng=%d\n",ng);
- fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m);
+ fprintf(ficgp,"# jk=1 to 2^%d=%d\n",nqveff,m);
for(jk=1; jk <=m; jk++) {
fprintf(ficgp,"# jk=%d\n",jk);
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
double *agemingood, *agemaxgood; /* Currently identical for all covariates */
- /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */
+ /* modcovmax=2*nqveff;/\* Max number of modalities. We suppose */
/* a covariate has 2 modalities, should be equal to ncovcombmax *\/ */
sumnewp = vector(1,ncovcombmax);
/************** Forecasting ******************/
-void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
+void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int nqveff){
/* proj1, year, month, day of starting projection
agemin, agemax range of age
dateprev1 dateprev2 range of dates during which prevalence is computed
printf("Computing forecasting: result on file '%s', please wait... \n", fileresf);
fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf);
- if (cptcoveff==0) ncodemax[cptcoveff]=1;
+ if (nqveff==0) ncodemax[nqveff]=1;
stepsize=(int) (stepm+YEARM-1)/YEARM;
if(jprojmean==0) jprojmean=1;
if(mprojmean==0) jprojmean=1;
- i1=cptcoveff;
+ i1=nqveff;
if (cptcovn < 1){i1=1;}
fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
/* if (h==(int)(YEARM*yearp)){ */
for(cptcov=1, k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+ for(cptcod=1;cptcod<=ncodemax[nqveff];cptcod++){
k=k+1;
fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficresf," yearproj age");
for (h=0; h<=nhstepm; h++){
if (h*hstepm/YEARM*stepm ==yearp) {
fprintf(ficresf,"\n");
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
}
}
/* /\************** Back Forecasting ******************\/ */
-/* void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ */
+/* void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int nqveff){ */
/* /\* back1, year, month, day of starting backection */
/* agemin, agemax range of age */
/* dateprev1 dateprev2 range of dates during which prevalence is computed */
/* printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
/* fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
-/* if (cptcoveff==0) ncodemax[cptcoveff]=1; */
+/* if (nqveff==0) ncodemax[nqveff]=1; */
/* /\* if (mobilav!=0) { *\/ */
/* /\* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
/* if(jprojmean==0) jprojmean=1; */
/* if(mprojmean==0) jprojmean=1; */
-/* i1=cptcoveff; */
+/* i1=nqveff; */
/* if (cptcovn < 1){i1=1;} */
/* fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); */
/* /\* if (h==(int)(YEARM*yearp)){ *\/ */
/* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
-/* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
+/* for(cptcod=1;cptcod<=ncodemax[nqveff];cptcod++){ */
/* k=k+1; */
/* fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#"); */
-/* for(j=1;j<=cptcoveff;j++) { */
+/* for(j=1;j<=nqveff;j++) { */
/* fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
/* } */
/* fprintf(ficresfb," yearbproj age"); */
/* for (h=0; h<=nhstepm; h++){ */
/* if (h*hstepm/YEARM*stepm ==yearp) { */
/* fprintf(ficresfb,"\n"); */
-/* for(j=1;j<=cptcoveff;j++) */
+/* for(j=1;j<=nqveff;j++) */
/* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
/* fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm); */
/* } */
printf("Computing forecasting: result on file '%s' \n", filerespop);
fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop);
- if (cptcoveff==0) ncodemax[cptcoveff]=1;
+ if (nqveff==0) ncodemax[nqveff]=1;
/* if (mobilav!=0) { */
/* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
}
for(cptcov=1,k=0;cptcov<=i2;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+ for(cptcod=1;cptcod<=ncodemax[nqveff];cptcod++){
k=k+1;
fprintf(ficrespop,"\n#******");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficrespop,"******\n");
while (*p1++ == *p2++);
}
-int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:
- * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age
- * - nagesqr = 1 if age*age in the model, otherwise 0.
- * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age
- * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age
- * - cptcovage number of covariates with age*products =2
- * - cptcovs number of simple covariates
- * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
- * which is a new column after the 9 (ncovcol) variables.
- * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual
- * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
- * Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
- * - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
- */
+int decodemodel ( char model[], int lastobs)
+ /**< This routine decode the model and returns:
+ * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age
+ * - nagesqr = 1 if age*age in the model, otherwise 0.
+ * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age
+ * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age
+ * - cptcovage number of covariates with age*products =2
+ * - cptcovs number of simple covariates
+ * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
+ * which is a new column after the 9 (ncovcol) variables.
+ * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual
+ * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
+ * Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
+ * - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
+ */
{
int i, j, k, ks;
int j1, k1, k2;
if ((strpt=strstr(model,"age*age")) !=0){
printf(" strpt=%s, model=%s\n",strpt, model);
if(strpt != model){
- printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model);
- fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model); fflush(ficlog);
- return 1;
- }
+ return 1;
+ }
nagesqr=1;
if (strstr(model,"+age*age") !=0)
- substrchaine(modelsav, model, "+age*age");
+ substrchaine(modelsav, model, "+age*age");
else if (strstr(model,"age*age+") !=0)
- substrchaine(modelsav, model, "age*age+");
+ substrchaine(modelsav, model, "age*age+");
else
- substrchaine(modelsav, model, "age*age");
+ substrchaine(modelsav, model, "age*age");
}else
nagesqr=0;
if (strlen(modelsav) >1){
j=nbocc(modelsav,'+'); /**< j=Number of '+' */
j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */
- cptcovs=j+1-j1; /**< Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =2 */
+ cptcovs=j+1-j1; /**< Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =5-3=2 */
cptcovt= j+1; /* Number of total covariates in the model, not including
- * cst, age and age*age
- * V1+V1*age+ V3 + V3*V4+age*age=> 4*/
- /* including age products which are counted in cptcovage.
- * but the covariates which are products must be treated
- * separately: ncovn=4- 2=2 (V1+V3). */
+ * cst, age and age*age
+ * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/
+ /* including age products which are counted in cptcovage.
+ * but the covariates which are products must be treated
+ * separately: ncovn=4- 2=2 (V1+V3). */
cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */
cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */
* k= 1 2 3 4 5 6 7 8
* cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
* covar[k,i], value of kth covariate if not including age for individual i:
- * covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8)
- * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[4]=3 Tvar[8]=8
+ * covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)
+ * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[2]=1 Tvar[4]=3 Tvar[8]=8
* if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and
* Tage[++cptcovage]=k
* if products, new covar are created after ncovcol with k1
If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
/* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
- printf("cptcovprod=%d ", cptcovprod);
- fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
-
- scanf("%d ",i);*/
-
+ printf("cptcovprod=%d ", cptcovprod);
+ fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
+
+ scanf("%d ",i);*/
+/* Dispatching in quantitative and time varying covariates */
+
+ for(k=1, ncoveff=0, nqveff=0, ntveff=0, nqtveff=0;k<=cptcovn; k++){ /* or cptocvt */
+ if (Tvar[k] <=ncovcol){
+ ncoveff++;
+ }else if( Tvar[k] <=ncovcol+nqv){
+ nqveff++;
+ }else if( Tvar[k] <=ncovcol+nqv+ntv){
+ ntveff++;
+ }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
+ nqtveff++;
+ }else
+ printf("Error in effective covariates \n");
+ }
return (0); /* with covar[new additional covariate if product] and Tage if age */
/*endread:*/
- printf("Exiting decodemodel: ");
- return (1);
+ printf("Exiting decodemodel: ");
+ return (1);
}
int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
agebase=ageminpar;
agelim=agemaxpar;
- i1=pow(2,cptcoveff);
+ i1=pow(2,ncoveff);
if (cptcovn < 1){i1=1;}
for(k=1; k<=i1;k++){
fprintf(ficrespl,"#******");
printf("#******");
fprintf(ficlog,"#******");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficrespl,"#Age ");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i);
/* for (age=agebase; age<=agebase; age++){ */
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
fprintf(ficrespl,"%.0f ",age );
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
tot=0.;
for(i=1; i<=nlstate;i++){
agelim=agemaxpar;
- i1=pow(2,cptcoveff);
+ i1=pow(2,nqveff);
if (cptcovn < 1){i1=1;}
for(k=1; k<=i1;k++){
fprintf(ficresplb,"#******");
printf("#******");
fprintf(ficlog,"#******");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficresplb,"#Age ");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
for(i=1; i<=nlstate;i++) fprintf(ficresplb," %d-%d ",i,i);
bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
}
fprintf(ficresplb,"%.0f ",age );
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
tot=0.;
for(i=1; i<=nlstate;i++){
/* hstepm=1; aff par mois*/
pstamp(ficrespij);
fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
- i1= pow(2,cptcoveff);
+ i1= pow(2,nqveff);
/* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
/* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
/* k=k+1; */
- for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ for (k=1; k <= (int) pow(2,nqveff); k++){
fprintf(ficrespij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrespij,"******\n");
/* hstepm=1; aff par mois*/
pstamp(ficrespijb);
fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
- i1= pow(2,cptcoveff);
+ i1= pow(2,nqveff);
/* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
/* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
/* k=k+1; */
- for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ for (k=1; k <= (int) pow(2,nqveff); k++){
fprintf(ficrespijb,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrespijb,"******\n");
if(invalidvarcomb[k]){
covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */
coqvar=matrix(1,nqv,1,n); /**< used in readdata */
cotvar=ma3x(1,maxwav,1,ntv,1,n); /**< used in readdata */
- cotqvar=ma3x(1,maxwav,1,ntqv,1,n); /**< used in readdata */
+ cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< used in readdata */
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
/* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
v1+v2*age+v2*v3 makes cptcovn = 3
printf("\n");
/*--------- results files --------------*/
- fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
+ fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, weightopt,model);
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
/*if((stepm == 1) && (strcmp(model,".")==0)){*/
if(prevfcast==1){
/* if(stepm ==1){*/
- prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
+ prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, nqveff);
}
if(backcast==1){
ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
/* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
- bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
+ bage, fage, firstpass, lastpass, anback2, p, nqveff); */
free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
- for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ for (k=1; k <= (int) pow(2,nqveff); k++){
fprintf(ficreseij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficreseij,"******\n");
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
- for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ for (k=1; k <= (int) pow(2,nqveff); k++){
fprintf(ficrest,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrest,"******\n");
fprintf(ficresstdeij,"\n#****** ");
fprintf(ficrescveij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++) {
+ for(j=1;j<=nqveff;j++) {
fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficrescveij,"******\n");
fprintf(ficresvij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficresvij,"******\n");
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
- for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ for (k=1; k <= (int) pow(2,nqveff); k++){
fprintf(ficresvpl,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
+ for(j=1;j<=nqveff;j++)
fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficresvpl,"******\n");
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
- free_ma3x(cotqvar,1,maxwav,1,ntqv,1,n);
+ free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
free_ma3x(cotvar,1,maxwav,1,ntv,1,n);
free_matrix(coqvar,1,maxwav,1,n);
free_matrix(covar,0,NCOVMAX,1,n);
git://henry.ined.fr / .git / commitdiff