/* $Id$
$State$
$Log$
+ Revision 1.222 2016/02/17 08:14:50 brouard
+ Summary: Probably last 0.98 stable version 0.98r6
+
Revision 1.221 2016/02/15 23:35:36 brouard
Summary: minor bug
/* #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 POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
+int nqv=0, ntv=0, nqtv=0; /* Total number of quantitative variables, time variable (dummy), quantitative and time variable */
int popbased=0;
int *wav; /* Number of waves for this individuual 0 is possible */
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 ***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;
/*double t34;*/
int i,j, nc, ii, jj;
- for(i=1; i<= nlstate; i++){
- for(j=1; j<i;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += param[i][j][nc]*cov[nc];*/
- lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
- /* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
- }
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- /* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
- }
- for(j=i+1; j<=nlstate+ndeath;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
- lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
- /* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
- }
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- }
- }
+ for(i=1; i<= nlstate; i++){
+ for(j=1; j<i;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += param[i][j][nc]*cov[nc];*/
+ lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+ /* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ /* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ for(j=i+1; j<=nlstate+ndeath;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+ lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+ /* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ }
+ }
- for(i=1; i<= nlstate; i++){
- s1=0;
- for(j=1; j<i; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- for(j=i+1; j<=nlstate+ndeath; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
- ps[i][i]=1./(s1+1.);
- /* Computing other pijs */
- for(j=1; j<i; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- for(j=i+1; j<=nlstate+ndeath; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
- } /* end i */
+ for(i=1; i<= nlstate; i++){
+ s1=0;
+ for(j=1; j<i; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ for(j=i+1; j<=nlstate+ndeath; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+ ps[i][i]=1./(s1+1.);
+ /* Computing other pijs */
+ for(j=1; j<i; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ for(j=i+1; j<=nlstate+ndeath; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+ } /* end i */
- for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
- for(jj=1; jj<= nlstate+ndeath; jj++){
- ps[ii][jj]=0;
- ps[ii][ii]=1;
- }
- }
+ for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ ps[ii][jj]=0;
+ ps[ii][ii]=1;
+ }
+ }
- /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
- /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
- /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
- /* } */
- /* printf("\n "); */
- /* } */
- /* printf("\n ");printf("%lf ",cov[2]);*/
- /*
- for(i=1; i<= npar; i++) printf("%f ",x[i]);
+ /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+ /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
+ /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+ /* } */
+ /* printf("\n "); */
+ /* } */
+ /* printf("\n ");printf("%lf ",cov[2]);*/
+ /*
+ for(i=1; i<= npar; i++) printf("%f ",x[i]);
goto end;*/
- return ps;
+ return ps;
}
/*************** backward transition probabilities ***************/
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;
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 variying (age dependent) and stores them in cov[]
+ 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[2+nagesqr+k]=covar[Tvar[k]][i];
+ 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]; */
+ }
+
/* 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);
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
cov[2]=agexact;
if(nagesqr==1)
- cov[3]= agexact*agexact;
+ 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 */
}
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 */
- mw[++mi][i]=m;
+ mw[++mi][i]=m;
}
if(m >=lastpass){
- if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
- if(firsthree == 0){
- printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
- 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);
- mw[++mi][i]=m;
- }
- if(s[m][i]==-2){ /* Vital status is really unknown */
- nbwarn++;
- if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */
- printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
- fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
- }
- break;
- }
- break;
+ 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);
+ 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);
+ mw[++mi][i]=m;
+ }
+ if(s[m][i]==-2){ /* Vital status is really unknown */
+ nbwarn++;
+ if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */
+ printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+ fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+ }
+ break;
+ }
+ break;
}
else
- m++;
+ m++;
}/* end while */
/* After last pass */
if (s[m][i] > nlstate){ /* In a death state */
mi++; /* Death is another wave */
/* if(mi==0) never been interviewed correctly before death */
- /* Only death is a correct wave */
+ /* Only death is a correct wave */
mw[mi][i]=m;
}else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */
/* m++; */
/* 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;
- }
- 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 );
+ 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;
+ }
+ 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){
- printf("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.\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 );
- firstfour=1;
- }
- 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 );
+ /* death is known but not confirmed by death status at any wave */
+ if(firstfour==0){
+ printf("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.\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 );
+ firstfour=1;
+ }
+ 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;
if(mi==0){
nbwarn++;
if(first==0){
- printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
- first=1;
+ printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
+ first=1;
}
if(first==1){
- fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
+ fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
}
} /* end mi==0 */
} /* End individuals */
/* wav and mw are no more changed */
-
+
for(i=1; i<=imx; i++){
for(mi=1; mi<wav[i];mi++){
if (stepm <=0)
- dh[mi][i]=1;
+ dh[mi][i]=1;
else{
- if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
- if (agedc[i] < 2*AGESUP) {
- j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);
- if(j==0) j=1; /* Survives at least one month after exam */
- else if(j<0){
- nberr++;
- printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- j=1; /* Temporary Dangerous patch */
- printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
- fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
- }
- k=k+1;
- if (j >= jmax){
- jmax=j;
- ijmax=i;
- }
- if (j <= jmin){
- jmin=j;
- ijmin=i;
- }
- sum=sum+j;
- /*if (j<0) printf("j=%d num=%d \n",j,i);*/
- /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
- }
- }
- else{
- j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
+ if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
+ if (agedc[i] < 2*AGESUP) {
+ j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);
+ if(j==0) j=1; /* Survives at least one month after exam */
+ else if(j<0){
+ nberr++;
+ printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ j=1; /* Temporary Dangerous patch */
+ printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
+ fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
+ }
+ k=k+1;
+ if (j >= jmax){
+ jmax=j;
+ ijmax=i;
+ }
+ if (j <= jmin){
+ jmin=j;
+ ijmin=i;
+ }
+ sum=sum+j;
+ /*if (j<0) printf("j=%d num=%d \n",j,i);*/
+ /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
+ }
+ }
+ else{
+ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
-
- k=k+1;
- if (j >= jmax) {
- jmax=j;
- ijmax=i;
- }
- else if (j <= jmin){
- jmin=j;
- ijmin=i;
- }
- /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
- /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
- if(j<0){
- nberr++;
- printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- }
- sum=sum+j;
- }
- jk= j/stepm;
- jl= j -jk*stepm;
- ju= j -(jk+1)*stepm;
- if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
- if(jl==0){
- dh[mi][i]=jk;
- bh[mi][i]=0;
- }else{ /* We want a negative bias in order to only have interpolation ie
- * to avoid the price of an extra matrix product in likelihood */
- dh[mi][i]=jk+1;
- bh[mi][i]=ju;
- }
- }else{
- if(jl <= -ju){
- dh[mi][i]=jk;
- bh[mi][i]=jl; /* bias is positive if real duration
- * is higher than the multiple of stepm and negative otherwise.
- */
- }
- else{
- dh[mi][i]=jk+1;
- bh[mi][i]=ju;
- }
- if(dh[mi][i]==0){
- dh[mi][i]=1; /* At least one step */
- bh[mi][i]=ju; /* At least one step */
- /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
- }
- } /* end if mle */
+
+ k=k+1;
+ if (j >= jmax) {
+ jmax=j;
+ ijmax=i;
+ }
+ else if (j <= jmin){
+ jmin=j;
+ ijmin=i;
+ }
+ /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
+ /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
+ if(j<0){
+ nberr++;
+ printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ }
+ sum=sum+j;
+ }
+ jk= j/stepm;
+ jl= j -jk*stepm;
+ ju= j -(jk+1)*stepm;
+ if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
+ if(jl==0){
+ dh[mi][i]=jk;
+ bh[mi][i]=0;
+ }else{ /* We want a negative bias in order to only have interpolation ie
+ * to avoid the price of an extra matrix product in likelihood */
+ dh[mi][i]=jk+1;
+ bh[mi][i]=ju;
+ }
+ }else{
+ if(jl <= -ju){
+ dh[mi][i]=jk;
+ bh[mi][i]=jl; /* bias is positive if real duration
+ * is higher than the multiple of stepm and negative otherwise.
+ */
+ }
+ else{
+ dh[mi][i]=jk+1;
+ bh[mi][i]=ju;
+ }
+ if(dh[mi][i]==0){
+ dh[mi][i]=1; /* At least one step */
+ bh[mi][i]=ju; /* At least one step */
+ /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
+ }
+ } /* end if mle */
}
} /* end wave */
}
}
/******************* Gnuplot file **************/
- void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
char dirfileres[132],optfileres[132];
- char gplotcondition[132];
+ char gplotcondition[132];
int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
int lv=0, vlv=0, kl=0;
int ng=0;
int vpopbased;
- int ioffset; /* variable offset for columns */
+ int ioffset; /* variable offset for columns */
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
/* printf("Problem with file %s",optionfilegnuplot); */
/*#ifdef windows */
fprintf(ficgp,"cd \"%s\" \n",pathc);
- /*#endif */
+ /*#endif */
m=pow(2,cptcoveff);
/* Contribution to likelihood */
/* Plot the probability implied in the likelihood */
- fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
- fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
- /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
- fprintf(ficgp,"\nset ter pngcairo size 640, 480");
+ fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
+ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
+ /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
+ fprintf(ficgp,"\nset ter pngcairo size 640, 480");
/* nice for mle=4 plot by number of matrix products.
replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */
/* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */
- /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
- fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
- fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
- fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
- fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
- for (i=1; i<= nlstate ; i ++) {
- fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
- fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk));
- fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
- for (j=2; j<= nlstate+ndeath ; j ++) {
- fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
- }
- fprintf(ficgp,";\nset out; unset ylabel;\n");
+ /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
+ fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
+ fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
+ fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
+ fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
+ for (i=1; i<= nlstate ; i ++) {
+ fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
+ fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk));
+ fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
+ for (j=2; j<= nlstate+ndeath ; j ++) {
+ fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
}
- /* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */
- /* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
- /* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
- fprintf(ficgp,"\nset out;unset log\n");
- /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
+ fprintf(ficgp,";\nset out; unset ylabel;\n");
+ }
+ /* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */
+ /* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
+ /* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
+ fprintf(ficgp,"\nset out;unset log\n");
+ /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
strcpy(dirfileres,optionfilefiname);
strcpy(optfileres,"vpl");
- /* 1eme*/
+ /* 1eme*/
for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
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 */
- /* 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 */
- vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
- /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ lv= decodtabm(k1,k,cptcoveff); /* 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 */
+ vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
+ /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
- fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
- fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
- fprintf(ficgp,"set xlabel \"Age\" \n\
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
+ fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \n\
set ylabel \"Probability\" \n \
set ter svg size 640, 480\n \
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
- for (i=1; i<= nlstate ; i ++) {
- if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
- for (i=1; i<= nlstate ; i ++) {
- if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
- for (i=1; i<= nlstate ; i ++) {
- if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
- 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 */
- 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 */
- /* 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 */
- vlv= nbcode[Tvaraff[k]][lv];
- kl++;
- /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
- /*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){
- 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{
- fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
- kl++;
- }
- } /* end covariate */
- }
- fprintf(ficgp,"\nset out \n");
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+ 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){
+ 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 */
+ /* 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 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ kl++;
+ /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+ /*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){
+ 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{
+ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+ kl++;
+ }
+ } /* end covariate */
+ } /* end if no covariate */
+ } /* end if backcast */
+ fprintf(ficgp,"\nset out \n");
} /* k1 */
} /* cpt */
/*2 eme*/
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 */
- /* 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 */
- vlv= nbcode[Tvaraff[k]][lv];
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
- }
- fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
+ 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 */
+ /* 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 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
- fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
- for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
- if(vpopbased==0)
- fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
- else
- fprintf(ficgp,"\nreplot ");
- for (i=1; i<= nlstate+1 ; i ++) {
- k=2*i;
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
- else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- fprintf(ficgp,"\" t\"\" w l lt 0,");
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
- else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
- } /* state */
- } /* vpopbased */
- fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+ fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+ for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+ if(vpopbased==0)
+ fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+ else
+ fprintf(ficgp,"\nreplot ");
+ for (i=1; i<= nlstate+1 ; i ++) {
+ k=2*i;
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+ else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"\" w l lt 0,");
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+ else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+ } /* state */
+ } /* vpopbased */
+ fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
} /* k1 */
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 */
- /* 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 */
- vlv= nbcode[Tvaraff[k]][lv];
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ lv= decodtabm(k1,k,cptcoveff); /* 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 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
/* k=2+nlstate*(2*cpt-2); */
k=2+(nlstate+1)*(cpt-1);
fprintf(ficgp,"set ter svg size 640, 480\n\
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
- for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
- fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
- fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
- for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
- fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+ fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
*/
for (i=1; i< nlstate ; i ++) {
- fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
- /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
+ /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
}
fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
}
}
- /* 4eme */
+ /* 4eme */
/* Survival functions (period) from state i in state j by initial state i */
for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
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 */
- /* 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 */
- vlv= nbcode[Tvaraff[k]][lv];
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ lv= decodtabm(k1,k,cptcoveff); /* 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 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
plot [%.f:%.f] ", ageminpar, agemaxpar);
k=3;
for (i=1; i<= nlstate ; i ++){
- if(i==1){
- fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
- }else{
- fprintf(ficgp,", '' ");
- }
- l=(nlstate+ndeath)*(i-1)+1;
- fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
- for (j=2; j<= nlstate+ndeath ; j ++)
- fprintf(ficgp,"+$%d",k+l+j-1);
- fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+ if(i==1){
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+ }else{
+ fprintf(ficgp,", '' ");
+ }
+ l=(nlstate+ndeath)*(i-1)+1;
+ fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+ for (j=2; j<= nlstate+ndeath ; j ++)
+ fprintf(ficgp,"+$%d",k+l+j-1);
+ fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
} /* nlstate */
fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
/* Survival functions (period) from state i in state j by final state j */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
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 */
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
/* CV preval stable (period) for each covariate */
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 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 */
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
-
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
+
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n\
-unset log y\n\
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
k=3; /* Offset */
for (i=1; i<= nlstate ; i ++){
fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
} /* end covariate */
-
-
+
+
/* 7eme */
if(backcast == 1){
/* CV back preval stable (period) for each covariate */
lv= decodtabm(k1,k,cptcoveff); /* 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 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
vlv= nbcode[Tvaraff[k]][lv];
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
}
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
} /* end covariate */
} /* End if backcast */
- /* 8eme */
+ /* 8eme */
if(prevfcast==1){
/* Projection from cross-sectional to stable (period) for each covariate */
}
fprintf(ficgp,"\n#\n");
if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
}
fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
-set ter svg size 640, 480\n \
-unset log y\n \
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
ioffset=4; /* Age is in 4 */
}else{
ioffset=6; /* Age is in 6 */
- /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
- /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
}
fprintf(ficgp," u %d:(",ioffset);
kl=0;
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
}else{
- fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
}
} /* end if covariate */
} /* nlstate */
fprintf(ficgp,"\nset out\n");
- } /* end cpt state*/
- } /* end covariate */
- } /* End if prevfcast */
+ } /* end cpt state*/
+ } /* end covariate */
+ } /* End if prevfcast */
- /* proba elementaires */
- fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+ /* proba elementaires */
+ fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
for(i=1,jk=1; i <=nlstate; i++){
fprintf(ficgp,"# initial state %d\n",i);
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
- fprintf(ficgp,"# current state %d\n",k);
- for(j=1; j <=ncovmodel; j++){
- fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
- jk++;
- }
- fprintf(ficgp,"\n");
+ fprintf(ficgp,"# current state %d\n",k);
+ for(j=1; j <=ncovmodel; j++){
+ fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
+ jk++;
+ }
+ fprintf(ficgp,"\n");
}
}
- }
+ }
fprintf(ficgp,"##############\n#\n");
-
+
/*goto avoid;*/
fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n");
fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
fprintf(ficgp,"#\n");
- for(ng=1; ng<=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);
- 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);
- fprintf(ficgp,"\nset ter svg size 640, 480 ");
- if (ng==1){
- fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
- fprintf(ficgp,"\nunset log y");
- }else if (ng==2){
- fprintf(ficgp,"\nset ylabel \"Probability\"\n");
- fprintf(ficgp,"\nset log y");
- }else if (ng==3){
- fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
- fprintf(ficgp,"\nset log y");
- }else
- fprintf(ficgp,"\nunset title ");
- fprintf(ficgp,"\nplot [%.f:%.f] ",ageminpar,agemaxpar);
- i=1;
- for(k2=1; k2<=nlstate; k2++) {
- k3=i;
- for(k=1; k<=(nlstate+ndeath); k++) {
- if (k != k2){
- switch( ng) {
- case 1:
- if(nagesqr==0)
- fprintf(ficgp," p%d+p%d*x",i,i+1);
- else /* nagesqr =1 */
- fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
- break;
- case 2: /* ng=2 */
- if(nagesqr==0)
- fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
- else /* nagesqr =1 */
- fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
- break;
- case 3:
- if(nagesqr==0)
- fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
- else /* nagesqr =1 */
- fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
- break;
- }
- ij=1;/* To be checked else nbcode[0][0] wrong */
- for(j=3; j <=ncovmodel-nagesqr; j++) {
- /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
- if(ij <=cptcovage) { /* Bug valgrind */
- if((j-2)==Tage[ij]) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
- /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
- ij++;
- }
- }
- else
- fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
- }
- }else{
- i=i-ncovmodel;
- if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
- fprintf(ficgp," (1.");
- }
+ 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);
+ 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);
+ fprintf(ficgp,"\nset ter svg size 640, 480 ");
+ if (ng==1){
+ fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
+ fprintf(ficgp,"\nunset log y");
+ }else if (ng==2){
+ fprintf(ficgp,"\nset ylabel \"Probability\"\n");
+ fprintf(ficgp,"\nset log y");
+ }else if (ng==3){
+ fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
+ fprintf(ficgp,"\nset log y");
+ }else
+ fprintf(ficgp,"\nunset title ");
+ fprintf(ficgp,"\nplot [%.f:%.f] ",ageminpar,agemaxpar);
+ i=1;
+ for(k2=1; k2<=nlstate; k2++) {
+ k3=i;
+ for(k=1; k<=(nlstate+ndeath); k++) {
+ if (k != k2){
+ switch( ng) {
+ case 1:
+ if(nagesqr==0)
+ fprintf(ficgp," p%d+p%d*x",i,i+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+ break;
+ case 2: /* ng=2 */
+ if(nagesqr==0)
+ fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+ break;
+ case 3:
+ if(nagesqr==0)
+ fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
+ break;
+ }
+ ij=1;/* To be checked else nbcode[0][0] wrong */
+ for(j=3; j <=ncovmodel-nagesqr; j++) {
+ /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
+ if(ij <=cptcovage) { /* Bug valgrind */
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+ ij++;
+ }
+ }
+ else
+ fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ }
+ }else{
+ i=i-ncovmodel;
+ if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
+ fprintf(ficgp," (1.");
+ }
- if(ng != 1){
- fprintf(ficgp,")/(1");
+ if(ng != 1){
+ fprintf(ficgp,")/(1");
- for(k1=1; k1 <=nlstate; k1++){
- if(nagesqr==0)
- fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
- else /* nagesqr =1 */
- fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+ for(k1=1; k1 <=nlstate; k1++){
+ if(nagesqr==0)
+ fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
- ij=1;
- for(j=3; j <=ncovmodel-nagesqr; j++){
- if(ij <=cptcovage) { /* Bug valgrind */
- if((j-2)==Tage[ij]) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
- /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
- ij++;
- }
- }
- else
- fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
- }
- fprintf(ficgp,")");
- }
- fprintf(ficgp,")");
- if(ng ==2)
- fprintf(ficgp," t \"p%d%d\" ", k2,k);
- else /* ng= 3 */
- fprintf(ficgp," t \"i%d%d\" ", k2,k);
- }else{ /* end ng <> 1 */
- if( k !=k2) /* logit p11 is hard to draw */
- fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
- }
- if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
- fprintf(ficgp,",");
- if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
- fprintf(ficgp,",");
- i=i+ncovmodel;
- } /* end k */
- } /* end k2 */
- fprintf(ficgp,"\n set out\n");
- } /* end jk */
- } /* end ng */
- /* avoid: */
- fflush(ficgp);
+ ij=1;
+ for(j=3; j <=ncovmodel-nagesqr; j++){
+ if(ij <=cptcovage) { /* Bug valgrind */
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+ ij++;
+ }
+ }
+ else
+ fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ }
+ fprintf(ficgp,")");
+ }
+ fprintf(ficgp,")");
+ if(ng ==2)
+ fprintf(ficgp," t \"p%d%d\" ", k2,k);
+ else /* ng= 3 */
+ fprintf(ficgp," t \"i%d%d\" ", k2,k);
+ }else{ /* end ng <> 1 */
+ if( k !=k2) /* logit p11 is hard to draw */
+ fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
+ }
+ if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
+ fprintf(ficgp,",");
+ if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
+ fprintf(ficgp,",");
+ i=i+ncovmodel;
+ } /* end k */
+ } /* end k2 */
+ fprintf(ficgp,"\n set out\n");
+ } /* end jk */
+ } /* end ng */
+ /* avoid: */
+ fflush(ficgp);
} /* end gnuplot */
/*-------- data file ----------*/
FILE *fic;
char dummy[]=" ";
- int i=0, j=0, n=0;
+ int i=0, j=0, n=0, iv=0;
+ int lstra;
int linei, month, year,iout;
char line[MAXLINE], linetmp[MAXLINE];
char stra[MAXLINE], strb[MAXLINE];
char *stratrunc;
- int lstra;
+
if((fic=fopen(datafile,"r"))==NULL) {
}
trimbb(linetmp,line); /* Trims multiple blanks in line */
strcpy(line, linetmp);
-
-
+
+ /* Loops on waves */
for (j=maxwav;j>=1;j--){
+ for (iv=nqtv;iv>=1;iv--){ /* Loop on time varying quantitative variables */
+ cutv(stra, strb, line, ' ');
+ if(strb[0]=='.') { /* Missing value */
+ lval=-1;
+ }else{
+ errno=0;
+ /* what_kind_of_number(strb); */
+ dval=strtod(strb,&endptr);
+ /* if( strb[0]=='\0' || (*endptr != '\0')){ */
+ /* if(strb != endptr && *endptr == '\0') */
+ /* dval=dlval; */
+ /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav);
+ fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog);
+ return 1;
+ }
+ cotqvar[j][iv][i]=dval;
+ }
+ strcpy(line,stra);
+ }/* end loop ntqv */
+
+ for (iv=ntv;iv>=1;iv--){ /* Loop on time varying dummies */
+ cutv(stra, strb, line, ' ');
+ if(strb[0]=='.') { /* Missing value */
+ lval=-1;
+ }else{
+ errno=0;
+ lval=strtol(strb,&endptr,10);
+ /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav);
+ fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog);
+ return 1;
+ }
+ }
+ if(lval <-1 || lval >1){
+ printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
+ For example, for multinomial values like 1, 2 and 3,\n \
+ build V1=0 V2=0 for the reference value (1),\n \
+ V1=1 V2=0 for (2) \n \
+ and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
+ output of IMaCh is often meaningless.\n \
+ Exiting.\n",lval,linei, i,line,j);
+ fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
+ For example, for multinomial values like 1, 2 and 3,\n \
+ build V1=0 V2=0 for the reference value (1),\n \
+ V1=1 V2=0 for (2) \n \
+ and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
+ output of IMaCh is often meaningless.\n \
+ Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
+ return 1;
+ }
+ cotvar[j][iv][i]=(double)(lval);
+ strcpy(line,stra);
+ }/* end loop ntv */
+
+ /* Statuses at wave */
cutv(stra, strb, line, ' ');
- if(strb[0]=='.') { /* Missing status */
- lval=-1;
+ if(strb[0]=='.') { /* Missing value */
+ lval=-1;
}else{
- errno=0;
- lval=strtol(strb,&endptr,10);
- /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
- fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
- return 1;
- }
+ errno=0;
+ lval=strtol(strb,&endptr,10);
+ /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
+ fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+ return 1;
+ }
}
+
s[j][i]=lval;
-
+
+ /* Date of Interview */
strcpy(line,stra);
cutv(stra, strb,line,' ');
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
}
else if( (iout=sscanf(strb,"%s.",dummy)) != 0){
- month=99;
- year=9999;
+ month=99;
+ year=9999;
}else{
- printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
- fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
- return 1;
+ printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
+ fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
+ return 1;
}
anint[j][i]= (double) year;
mint[j][i]= (double)month;
strcpy(line,stra);
- } /* ENd Waves */
-
+ } /* End loop on waves */
+
+ /* Date of death */
cutv(stra, strb,line,' ');
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
}
year=9999;
}else{
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
- fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog);
- return 1;
+ fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog);
+ return 1;
}
andc[i]=(double) year;
moisdc[i]=(double) month;
strcpy(line,stra);
+ /* Date of birth */
cutv(stra, strb,line,' ');
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
}
}else{
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog);
- return 1;
+ return 1;
}
if (year==9999) {
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);
- return 1;
+ return 1;
}
annais[i]=(double)(year);
moisnais[i]=(double)(month);
strcpy(line,stra);
-
+
+ /* Sample weight */
cutv(stra, strb,line,' ');
errno=0;
dval=strtod(strb,&endptr);
}
weight[i]=dval;
strcpy(line,stra);
+
+ for (iv=nqv;iv>=1;iv--){ /* Loop on fixed quantitative variables */
+ cutv(stra, strb, line, ' ');
+ if(strb[0]=='.') { /* Missing value */
+ lval=-1;
+ }else{
+ errno=0;
+ /* what_kind_of_number(strb); */
+ dval=strtod(strb,&endptr);
+ /* if(strb != endptr && *endptr == '\0') */
+ /* dval=dlval; */
+ /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value (out of %d) constant for all waves. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqv, maxwav);
+ fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value (out of %d) constant for all waves. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqv, maxwav);fflush(ficlog);
+ return 1;
+ }
+ coqvar[iv][i]=dval;
+ }
+ strcpy(line,stra);
+ }/* end loop nqv */
+ /* Covariate values */
for (j=ncovcol;j>=1;j--){
cutv(stra, strb,line,' ');
- if(strb[0]=='.') { /* Missing status */
- lval=-1;
+ if(strb[0]=='.') { /* Missing covariate value */
+ lval=-1;
}else{
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);
- fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog);
- return 1;
- }
+ errno=0;
+ lval=strtol(strb,&endptr,10);
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);
+ fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog);
+ return 1;
+ }
}
if(lval <-1 || lval >1){
- printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+ printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
For example, for multinomial values like 1, 2 and 3,\n \
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
output of IMaCh is often meaningless.\n \
Exiting.\n",lval,linei, i,line,j);
- fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+ fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
For example, for multinomial values like 1, 2 and 3,\n \
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
output of IMaCh is often meaningless.\n \
Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
- return 1;
+ return 1;
}
covar[j][i]=(double)(lval);
strcpy(line,stra);
return (0);
/* endread: */
- printf("Exiting readdata: ");
- fclose(fic);
- return (1);
-
-
-
+ printf("Exiting readdata: ");
+ fclose(fic);
+ return (1);
}
+
void removespace(char *str) {
char *p1 = str, *p2 = str;
do
Tvar[k]=0;
cptcovage=0;
for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
- cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
- modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */
- if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
- /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
- /*scanf("%d",i);*/
- if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
- cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
- if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
- /* covar is not filled and then is empty */
- cptcovprod--;
- cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
- Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
- cptcovage++; /* Sums the number of covariates which include age as a product */
- Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
- /*printf("stre=%s ", stre);*/
- } else if (strcmp(strd,"age")==0) { /* or age*Vn */
- cptcovprod--;
- cutl(stre,strb,strc,'V');
- Tvar[k]=atoi(stre);
- cptcovage++;
- Tage[cptcovage]=k;
- } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/
- /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
- cptcovn++;
- cptcovprodnoage++;k1++;
- cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
- Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
- because this model-covariate is a construction we invent a new column
- ncovcol + k1
- If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
- Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
- cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
- Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */
- Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
- Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
- k2=k2+2;
- Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
- Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
- for (i=1; i<=lastobs;i++){
- /* Computes the new covariate which is a product of
- covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
- covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
- }
- } /* End age is not in the model */
- } /* End if model includes a product */
- else { /* no more sum */
- /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
- /* scanf("%d",i);*/
- cutl(strd,strc,strb,'V');
- ks++; /**< Number of simple covariates */
- cptcovn++;
- Tvar[k]=atoi(strd);
- }
- strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
- /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
- scanf("%d",i);*/
+ cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
+ modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */
+ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
+ /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
+ /*scanf("%d",i);*/
+ if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
+ cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
+ if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
+ /* covar is not filled and then is empty */
+ cptcovprod--;
+ cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
+ Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
+ cptcovage++; /* Sums the number of covariates which include age as a product */
+ Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
+ /*printf("stre=%s ", stre);*/
+ } else if (strcmp(strd,"age")==0) { /* or age*Vn */
+ cptcovprod--;
+ cutl(stre,strb,strc,'V');
+ Tvar[k]=atoi(stre);
+ cptcovage++;
+ Tage[cptcovage]=k;
+ } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/
+ /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+ cptcovn++;
+ cptcovprodnoage++;k1++;
+ cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+ Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
+ because this model-covariate is a construction we invent a new column
+ ncovcol + k1
+ If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
+ Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
+ cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
+ Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */
+ Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
+ Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
+ k2=k2+2;
+ Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
+ Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
+ for (i=1; i<=lastobs;i++){
+ /* Computes the new covariate which is a product of
+ covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+ covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+ }
+ } /* End age is not in the model */
+ } /* End if model includes a product */
+ else { /* no more sum */
+ /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
+ /* scanf("%d",i);*/
+ cutl(strd,strc,strb,'V');
+ ks++; /**< Number of simple covariates */
+ cptcovn++;
+ Tvar[k]=atoi(strd);
+ }
+ strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
+ /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
+ scanf("%d",i);*/
} /* end of loop + on total covariates */
} /* end if strlen(modelsave == 0) age*age might exist */
} /* end if strlen(model == 0) */
}else
break;
}
- if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
- &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
- if (num_filled != 8) {
- printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
+ if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
+ &ftol, &stepm, &ncovcol, &nqv, &ntv, &nqtv, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
+ if (num_filled != 11) {
+ printf("Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
printf("but line=%s\n",line);
}
- printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
+ printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt);
}
/* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
/*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
/* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
/* numlinepar=numlinepar+3; /\* In general *\/ */
/* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
- fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
- fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+ fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model);
+ fprintf(ficlog,"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=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model);
fflush(ficlog);
/* if(model[0]=='#'|| model[0]== '\0'){ */
if(model[0]=='#'){
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 */
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
/* Main decodemodel */
- if(decodemodel(model, lastobs) == 1)
+ if(decodemodel(model, lastobs) == 1) /* In order to get Tvar[k] V4+V3+V5 p Tvar[1]@3 = {4, 3, 5}*/
goto end;
if((double)(lastobs-imx)/(double)imx > 1.10){
and for any valid combination of covariates
and prints on file fileres'p'. */
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
- firstpass, lastpass, stepm, weightopt, model);
+ firstpass, lastpass, stepm, weightopt, model);
fprintf(fichtm,"\n");
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
ageexmed=vector(1,n);
agecens=vector(1,n);
dcwave=ivector(1,n);
-
+
for (i=1; i<=imx; i++){
dcwave[i]=-1;
for (m=firstpass; m<=lastpass; m++)
ungetc(c,ficpar);
fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
- fprintf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
- fprintf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
- fprintf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+ fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+ fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+ fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
/* day and month of proj2 are not used but only year anproj2.*/
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(cotvar,1,maxwav,1,ntv,1,n);
+ free_matrix(coqvar,1,maxwav,1,n);
free_matrix(covar,0,NCOVMAX,1,n);
free_matrix(matcov,1,npar,1,npar);
free_matrix(hess,1,npar,1,npar);
git://henry.ined.fr / .git / commitdiff