--- imach/src/imach.c 2018/02/27 22:50:02 1.282 +++ imach/src/imach.c 2019/05/09 13:39:37 1.290 @@ -1,6 +1,32 @@ -/* $Id: imach.c,v 1.282 2018/02/27 22:50:02 brouard Exp $ +/* $Id: imach.c,v 1.290 2019/05/09 13:39:37 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.290 2019/05/09 13:39:37 brouard + Summary: 0.99r18 unlimited number of individuals + + The number n which was limited to 20,000 cases is now unlimited, from firstobs to lastobs. If the number is too for the virtual memory, probably an error will occur. + + Revision 1.289 2018/12/13 09:16:26 brouard + Summary: Bug for young ages (<-30) will be in r17 + + Revision 1.288 2018/05/02 20:58:27 brouard + Summary: Some bugs fixed + + Revision 1.287 2018/05/01 17:57:25 brouard + Summary: Bug fixed by providing frequencies only for non missing covariates + + Revision 1.286 2018/04/27 14:27:04 brouard + Summary: some minor bugs + + Revision 1.285 2018/04/21 21:02:16 brouard + Summary: Some bugs fixed, valgrind tested + + Revision 1.284 2018/04/20 05:22:13 brouard + Summary: Computing mean and stdeviation of fixed quantitative variables + + Revision 1.283 2018/04/19 14:49:16 brouard + Summary: Some minor bugs fixed + Revision 1.282 2018/02/27 22:50:02 brouard *** empty log message *** @@ -1028,14 +1054,15 @@ typedef struct { #define NINTERVMAX 8 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ -#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ +/* #define NCOVMAX 20 */ /**< Maximum number of covariates, including generated covariates V1*V2 */ #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 /*#define decodtabm(h,k,cptcoveff)= (h <= (1<> (k-1)) & 1) +1 : -1)*/ #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 -#define MAXN 20000 +/*#define MAXN 20000 */ /* Should by replaced by nobs, real number of observations and unlimited */ #define YEARM 12. /**< Number of months per year */ /* #define AGESUP 130 */ -#define AGESUP 150 +/* #define AGESUP 150 */ +#define AGESUP 200 #define AGEINF 0 #define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */ #define AGEBASE 40 @@ -1051,12 +1078,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.282 2018/02/27 22:50:02 brouard Exp $ */ +/* $Id: imach.c,v 1.290 2019/05/09 13:39:37 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -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: 1.282 $ $Date: 2018/02/27 22:50:02 $"; +char copyright[]="April 2018,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: 1.290 $ $Date: 2019/05/09 13:39:37 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1079,6 +1106,7 @@ int nqfveff=0; /**< nqfveff Number of Qu 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 nobs=10; /* Number of observations in the data lastobs-firstobs */ int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */ int npar=NPARMAX; int nlstate=2; /* Number of live states */ @@ -2562,6 +2590,7 @@ void powell(double p[], double **xi, int double **newm; double agefin, delaymax=200. ; /* 100 Max number of years to converge */ int ncvloop=0; + int first=0; min=vector(1,nlstate); max=vector(1,nlstate); @@ -2658,10 +2687,14 @@ void powell(double p[], double **xi, int free_vector(meandiff,1,nlstate); return prlim; } - } /* age loop */ + } /* agefin loop */ /* After some age loop it doesn't converge */ - printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ -Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + if(!first){ + first=1; + printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d). Others in log file only...\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); + } + fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); + /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ free_vector(min,1,nlstate); free_vector(max,1,nlstate); @@ -2727,7 +2760,8 @@ Earliest age to start was %d-%d=%d, ncvl /* Even if hstepm = 1, at least one multiplication by the unit matrix */ /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */ /* for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */ - for(agefin=age; agefin ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\ Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); @@ -4368,7 +4402,7 @@ void freqsummary(char fileres[], double double ***freq; /* Frequencies */ double *x, *y, a=0.,b=0.,r=1., sa=0., sb=0.; /* for regression, y=b+m*x and r is the correlation coefficient */ int no=0, linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb); - double *meanq; + double *meanq, *stdq, *idq; double **meanqt; double *pp, **prop, *posprop, *pospropt; double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; @@ -4381,6 +4415,8 @@ void freqsummary(char fileres[], double pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ + stdq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ + idq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ meanqt=matrix(1,lastpass,1,nqtveff); strcpy(fileresp,"P_"); strcat(fileresp,fileresu); @@ -4489,13 +4525,16 @@ Title=%s
Datafile=%s Firstpass=%d La posprop[i]=0; pospropt[i]=0; } - /* for (z1=1; z1<= nqfveff; z1++) { */ - /* meanq[z1]+=0.; */ + for (z1=1; z1<= nqfveff; z1++) { /* zeroing for each combination j1 as well as for the total */ + idq[z1]=0.; + meanq[z1]=0.; + stdq[z1]=0.; + } + /* for (z1=1; z1<= nqtveff; z1++) { */ /* for(m=1;m<=lastpass;m++){ */ - /* meanqt[m][z1]=0.; */ - /* } */ - /* } */ - + /* meanqt[m][z1]=0.; */ + /* } */ + /* } */ /* dateintsum=0; */ /* k2cpt=0; */ @@ -4505,9 +4544,6 @@ Title=%s
Datafile=%s Firstpass=%d La if(j !=0){ if(anyvaryingduminmodel==0){ /* If All fixed covariates */ if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - /* for (z1=1; z1<= nqfveff; z1++) { */ - /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ - /* } */ for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ /* if(Tvaraff[z1] ==-20){ */ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ @@ -4528,7 +4564,7 @@ Title=%s
Datafile=%s Firstpass=%d La }/* end j==0 */ if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ /* for(m=firstpass; m<=lastpass; m++){ */ - for(mi=1; miDatafile=%s Firstpass=%d La }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ } /* end j==0 */ /* bool =0 we keep that guy which corresponds to the combination of dummy values */ - if(bool==1){ + if(bool==1){ /*Selected */ /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] and mw[mi+1][iind]. dh depends on stepm. */ agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ @@ -4566,6 +4602,11 @@ Title=%s
Datafile=%s Firstpass=%d La if(s[m][iind]==-1) printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.)); freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ + for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */ + idq[z1]=idq[z1]+weight[iind]; + meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /* Computes mean of quantitative with selected filter */ + stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */ + } /* if((int)agev[m][iind] == 55) */ /* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ @@ -4581,6 +4622,11 @@ Title=%s
Datafile=%s Firstpass=%d La bool=1; }/* end bool 2 */ } /* end m */ + /* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */ + /* idq[z1]=idq[z1]+weight[iind]; */ + /* meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /\* Computes mean of quantitative with selected filter *\/ */ + /* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/ /\* Computes mean of quantitative with selected filter *\/ */ + /* } */ } /* end bool */ } /* end iind = 1 to imx */ /* prop[s][age] is feeded for any initial and valid live state as well as @@ -4618,6 +4664,27 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtmfr, "**********\n"); fprintf(ficlog, "**********\n"); } + /* + Printing means of quantitative variables if any + */ + for (z1=1; z1<= nqfveff; z1++) { + fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]); + fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]); + if(weightopt==1){ + printf(" Weighted mean and standard deviation of"); + fprintf(ficlog," Weighted mean and standard deviation of"); + fprintf(ficresphtmfr," Weighted mean and standard deviation of"); + } + printf(" fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + fprintf(ficlog," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + fprintf(ficresphtmfr," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)

\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + } + /* for (z1=1; z1<= nqtveff; z1++) { */ + /* for(m=1;m<=lastpass;m++){ */ + /* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f

\n", z1, m, meanqt[m][z1]); */ + /* } */ + /* } */ + fprintf(ficresphtm,""); if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ fprintf(ficresp, " Age"); @@ -4852,7 +4919,7 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficlog,"\n"); } } - } + } /* end of state i */ printf("#Freqsummary\n"); fprintf(ficlog,"\n"); for(s1=-1; s1 <=nlstate+ndeath; s1++){ @@ -4898,7 +4965,9 @@ Title=%s
Datafile=%s Firstpass=%d La fclose(ficresp); fclose(ficresphtm); fclose(ficresphtmfr); + free_vector(idq,1,nqfveff); free_vector(meanq,1,nqfveff); + free_vector(stdq,1,nqfveff); free_matrix(meanqt,1,lastpass,1,nqtveff); free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); @@ -5005,7 +5074,7 @@ void prevalence(double ***probs, double /*j=cptcoveff;*/ if (cptcovn<1) {j=1;ncodemax[1]=1;} - first=1; + first=0; for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ for (i=1; i<=nlstate; i++) for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) @@ -5063,12 +5132,11 @@ void prevalence(double ***probs, double if(posprop>=1.e-5){ probs[i][jk][j1]= prop[jk][i]/posprop; } else{ - if(first==1){ - first=0; + if(!first){ + first=1; printf("Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); - fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); }else{ - fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); + fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases.\n",jk,i,jk, j1,probs[i][jk][j1]); } } } @@ -5314,9 +5382,11 @@ void concatwav(int wav[], int **dh, int /* *cptcov=0; */ for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ + for (k=1; k <= maxncov; k++) + for(j=1; j<=2; j++) + nbcode[k][j]=0; /* Valgrind */ /* 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 (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ for (j=-1; (j < maxncov); j++) Ndum[j]=0; if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ @@ -5334,7 +5404,11 @@ void concatwav(int wav[], int **dh, int modmaxcovj=ij; else if (ij < modmincovj) modmincovj=ij; - if ((ij < -1) && (ij > NCOVMAX)){ + if (ij <0 || ij >1 ){ + printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); + fprintf(ficlog,"Information, currently IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); + } + if ((ij < -1) || (ij > NCOVMAX)){ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); exit(1); }else @@ -5380,12 +5454,18 @@ void concatwav(int wav[], int **dh, int /* nbcode[Tvar[j]][3]=2; */ /* To be continued (not working yet). */ ij=0; /* ij is similar to i but can jump over null modalities */ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + + /* for (i=modmincovj; i<=modmaxcovj; i++) { */ /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + /* Skipping the case of missing values by reducing nbcode to 0 and 1 and not -1, 0, 1 */ + /* model=V1+V2+V3, if V2=-1, 0 or 1, then nbcode[2][1]=0 and nbcode[2][2]=1 instead of + * nbcode[2][1]=-1, nbcode[2][2]=0 and nbcode[2][3]=1 */ + /*, could be restored in the future */ + for (i=0; i<=1; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ if (Ndum[i] == 0) { /* If nobody responded to this modality k */ break; } ij++; - nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/ + nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1 . Could be -1*/ cptcode = ij; /* New max modality for covar j */ } /* end of loop on modality i=-1 to 1 or more */ break; @@ -5401,21 +5481,7 @@ void concatwav(int wav[], int **dh, int break; } /* end switch */ } /* end dummy test */ - - /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ - /* /\*recode from 0 *\/ */ - /* k is a modality. If we have model=V1+V1*sex */ - /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ - /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ - /* } */ - /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ - /* if (ij > ncodemax[j]) { */ - /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ - /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ - /* break; */ - /* } */ - /* } /\* end of loop on modality k *\/ */ - } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ + } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; /* Look at fixed dummy (single or product) covariates to check empty modalities */ @@ -5799,6 +5865,7 @@ void concatwav(int wav[], int **dh, int double **dnewm,**doldm; double **dnewmp,**doldmp; int i, j, nhstepm, hstepm, h, nstepm ; + int first=0; int k; double *xp; double **gp, **gm; /**< for var eij */ @@ -5923,7 +5990,7 @@ void concatwav(int wav[], int **dh, int } /**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and * returns into prlim . - */ + */ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres); /* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */ @@ -5961,7 +6028,7 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres); if (popbased==1) { @@ -6147,7 +6214,7 @@ void concatwav(int wav[], int **dh, int int theta; pstamp(ficresvpl); - fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); + fprintf(ficresvpl,"# Standard deviation of period (forward stable) prevalences \n"); fprintf(ficresvpl,"# Age "); if(nresult >=1) fprintf(ficresvpl," Result# "); @@ -6176,20 +6243,20 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); - else - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ + /* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ + /* else */ + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); for(i=1;i<=nlstate;i++){ gp[i] = prlim[i][i]; mgp[theta][i] = prlim[i][i]; } for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); - else - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ + /* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ + /* else */ + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); for(i=1;i<=nlstate;i++){ gm[i] = prlim[i][i]; mgm[theta][i] = prlim[i][i]; @@ -6238,8 +6305,11 @@ void concatwav(int wav[], int **dh, int fprintf(ficresvpl,"%.0f ",age ); if(nresult >=1) fprintf(ficresvpl,"%d ",nres ); - for(i=1; i<=nlstate;i++) + for(i=1; i<=nlstate;i++){ fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); + /* for(j=1;j<=nlstate;j++) */ + /* fprintf(ficresvpl," %d %.5f ",j,prlim[j][i]); */ + } fprintf(ficresvpl,"\n"); free_vector(gp,1,nlstate); free_vector(gm,1,nlstate); @@ -6456,7 +6526,7 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

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

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

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. %s
    • \n",optionfilehtmcov,optionfilehtmcov); + fprintf(fichtm,"\n

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

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. File %s
    • \n",optionfilehtmcov,optionfilehtmcov); fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ and drawn. It helps understanding how is the covariance between two incidences.\ @@ -6749,10 +6819,10 @@ void printinghtml(char fileresu[], char - Estimated back transition probabilities over %d (stepm) months: %s
    \n ", stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_")); fprintf(fichtm,"\ - - Period (stable) prevalence in each health state: %s
    \n", + - Period (forward) prevalence in each health state: %s
    \n", subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); fprintf(fichtm,"\ - - Period (stable) back prevalence in each health state: %s
    \n", + - Backward prevalence in each health state: %s
    \n", subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_")); fprintf(fichtm,"\ - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ @@ -6865,24 +6935,24 @@ divided by h: hPij for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
    \n- Survival functions from state %d in each live state and total.\ Or probability to survive in various states (1 to %d) being in state %d at different ages. \ - %s_%d%d-%d.svg
    ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); + %s_%d-%d-%d.svg
    ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); } - /* Period (stable) prevalence in each health state */ + /* Period (forward stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
    \n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
    \ ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); } if(backcast==1){ - /* Period (stable) back prevalence in each health state */ + /* Backward prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
    \n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. %s_%d-%d-%d.svg
    \ ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); } } if(prevfcast==1){ - /* Projection of prevalence up to period (stable) prevalence in each health state */ + /* Projection of prevalence up to period (forward stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
    \n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
    \ + fprintf(fichtm,"
    \n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
    \ ", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); } } @@ -6935,13 +7005,13 @@ See page 'Matrix of variance-covariance %s
    \n", estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); fprintf(fichtm,"\ - - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
    \n", + - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the forward (period) prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
    \n", estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_")); fprintf(fichtm,"\ - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
    \n", estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_")); fprintf(fichtm,"\ - - Standard deviation of period (stable) prevalences: %s
    \n",\ + - Standard deviation of forward (period) prevalences: %s
    \n",\ subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); /* if(popforecast==1) fprintf(fichtm,"\n */ @@ -7077,7 +7147,7 @@ void printinggnuplot(char fileresu[], ch continue; /* We are interested in selected combination by the resultline */ /* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */ - fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); + fprintf(ficgp,"\n# 1st: Forward (stable period) prevalence with CI: 'VPL_' files and live state =%d ", cpt); strcpy(gplotlabel,"("); 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 */ @@ -7115,7 +7185,7 @@ void printinggnuplot(char fileresu[], ch 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==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); + fprintf(ficgp,"\" t\"Forward prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); @@ -7444,7 +7514,7 @@ set ter svg size 640, 480\nunset log y\n continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */ strcpy(gplotlabel,"("); - fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n#CV preval stable (forward): '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 */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -7488,14 +7558,14 @@ set ter svg size 640, 480\nunset log y\n /* 7eme */ if(backcast == 1){ - /* CV back preval stable (period) for each covariate */ + /* CV backward prevalence for each covariate */ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ if(m != 1 && TKresult[nres]!= k1) continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */ strcpy(gplotlabel,"("); - fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n#CV Backward stable prevalence: 'pijb' 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 */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -7543,7 +7613,7 @@ set ter svg size 640, 480\nunset log y\n /* 8eme */ if(prevfcast==1){ - /* Projection from cross-sectional to stable (period) for each covariate */ + /* Projection from cross-sectional to forward stable (period) prevalence for each covariate */ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ @@ -7551,7 +7621,7 @@ set ter svg size 640, 480\nunset log y\n continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ strcpy(gplotlabel,"("); - fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n#Projection of prevalence to forward stable prevalence (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 */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -7985,6 +8055,7 @@ set ter svg size 640, 480\nunset log y\n int modcovmax =1; int mobilavrange, mob; int iage=0; + int firstA1=0, firstA2=0; double sum=0., sumr=0.; double age; @@ -8082,6 +8153,7 @@ set ter svg size 640, 480\nunset log y\n } /* age */ /* Thus we have agemingood and agemaxgood as well as goodr for raw (preobs) */ /* but they will change */ + firstA1=0;firstA2=0; for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, filling up to the youngest */ sumnewm[cptcod]=0.; sumnewmr[cptcod]=0.; @@ -8114,11 +8186,19 @@ set ter svg size 640, 480\nunset log y\n sumr+=probs[(int)age][i][cptcod]; } if(fabs(sum - 1.) > 1.e-3) { /* bad */ - printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); + if(!firstA1){ + firstA1=1; + printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage); + } + fprintf(ficlog,"Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); } /* end bad */ /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ if(fabs(sumr - 1.) > 1.e-3) { /* bad */ - printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); + if(!firstA2){ + firstA2=1; + printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage); + } + fprintf(ficlog,"Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); } /* end bad */ }/* age */ @@ -8493,7 +8573,7 @@ set ter svg size 640, 480\nunset log y\n /* Variance of prevalence limit: varprlim */ void varprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **prlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){ - /*------- Variance of period (stable) prevalence------*/ + /*------- Variance of forward period (stable) prevalence------*/ char fileresvpl[FILENAMELENGTH]; FILE *ficresvpl; @@ -8504,11 +8584,11 @@ set ter svg size 640, 480\nunset log y\n strcpy(fileresvpl,"VPL_"); strcat(fileresvpl,fileresu); if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { - printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); + printf("Problem with variance of forward period (stable) prevalence resultfile: %s\n", fileresvpl); exit(0); } - printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); - fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); + printf("Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); + fprintf(ficlog, "Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ @@ -8545,8 +8625,8 @@ set ter svg size 640, 480\nunset log y\n } fclose(ficresvpl); - printf("done variance-covariance of period prevalence\n");fflush(stdout); - fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); + printf("done variance-covariance of forward period prevalence\n");fflush(stdout); + fprintf(ficlog,"done variance-covariance of forward period prevalence\n");fflush(ficlog); } /* Variance of back prevalence: varbprlim */ @@ -9689,7 +9769,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); - for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} + for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ Fixed[k]= 0; @@ -9939,11 +10019,12 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( /* Searching for doublons in the model */ for(k1=1; k1<= cptcovt;k1++){ for(k2=1; k2 maxwav){ + printf("Error (lastpass = %d) > (maxwav = %d)\n",lastpass, maxwav); + fprintf(ficlog,"Error (lastpass = %d) > (maxwav = %d)\n",lastpass, maxwav); + fflush(ficlog); + goto end; + } + 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); + fprintf(ficparo,"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); + fprintf(ficres,"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, 0, weightopt); + fprintf(ficlog,"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 */ @@ -10980,10 +11078,10 @@ int main(int argc, char *argv[]) /* If line starts with a # it is a comment */ if (line[0] == '#') { numlinepar++; - fputs(line,stdout); - fputs(line,ficparo); - fputs(line,ficres); - fputs(line,ficlog); + printf("%s",line); + fprintf(ficres,"%s",line); + fprintf(ficparo,"%s",line); + fprintf(ficlog,"%s",line); continue; }else break; @@ -11004,12 +11102,15 @@ int main(int argc, char *argv[]) } /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */ printf("model=1+age+%s\n",model);fflush(stdout); + fprintf(ficparo,"model=1+age+%s\n",model);fflush(stdout); + fprintf(ficres,"model=1+age+%s\n",model);fflush(stdout); + fprintf(ficlog,"model=1+age+%s\n",model);fflush(stdout); } /* 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 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); + /* 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]=='#'){ @@ -11037,10 +11138,10 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); - covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ - if(nqv>=1)coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */ - if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */ - if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/ + covar=matrix(0,NCOVMAX,firstobs,lastobs); /**< used in readdata */ + if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs); /**< Fixed quantitative covariate */ + if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs); /**< Time varying quantitative covariate */ + if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs); /**< Time varying covariate (dummy and quantitative)*/ 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 @@ -11243,16 +11344,25 @@ Please run with mle=-1 to get a correct /* Main data */ - n= lastobs; - num=lvector(1,n); - moisnais=vector(1,n); - annais=vector(1,n); - moisdc=vector(1,n); - andc=vector(1,n); - weight=vector(1,n); - agedc=vector(1,n); - cod=ivector(1,n); - for(i=1;i<=n;i++){ + nobs=lastobs-firstobs+1; /* was = lastobs;*/ + /* num=lvector(1,n); */ + /* moisnais=vector(1,n); */ + /* annais=vector(1,n); */ + /* moisdc=vector(1,n); */ + /* andc=vector(1,n); */ + /* weight=vector(1,n); */ + /* agedc=vector(1,n); */ + /* cod=ivector(1,n); */ + /* for(i=1;i<=n;i++){ */ + num=lvector(firstobs,lastobs); + moisnais=vector(firstobs,lastobs); + annais=vector(firstobs,lastobs); + moisdc=vector(firstobs,lastobs); + andc=vector(firstobs,lastobs); + weight=vector(firstobs,lastobs); + agedc=vector(firstobs,lastobs); + cod=ivector(firstobs,lastobs); + for(i=firstobs;i<=lastobs;i++){ num[i]=0; moisnais[i]=0; annais[i]=0; @@ -11262,9 +11372,9 @@ Please run with mle=-1 to get a correct cod[i]=0; weight[i]=1.0; /* Equal weights, 1 by default */ } - mint=matrix(1,maxwav,1,n); - anint=matrix(1,maxwav,1,n); - s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ + mint=matrix(1,maxwav,firstobs,lastobs); + anint=matrix(1,maxwav,firstobs,lastobs); + s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ tab=ivector(1,NCOVMAX); ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ @@ -11366,8 +11476,8 @@ Please run with mle=-1 to get a correct agegomp=(int)agemin; - free_vector(moisnais,1,n); - free_vector(annais,1,n); + free_vector(moisnais,firstobs,lastobs); + free_vector(annais,firstobs,lastobs); /* free_matrix(mint,1,maxwav,1,n); free_matrix(anint,1,maxwav,1,n);*/ /* free_vector(moisdc,1,n); */ @@ -11393,8 +11503,8 @@ Please run with mle=-1 to get a correct concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); /* Concatenates waves */ - free_vector(moisdc,1,n); - free_vector(andc,1,n); + free_vector(moisdc,firstobs,lastobs); + free_vector(andc,firstobs,lastobs); /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); @@ -11576,7 +11686,7 @@ Title=%s
    Datafile=%s Firstpass=%d La firstpass, lastpass, stepm, weightopt, model); fprintf(fichtm,"\n"); - fprintf(fichtm,"

    Parameter line 2

    • Tolerance for the convergence of the likelihood: ftol=%f \n
    • Interval for the elementary matrix (in month): stepm=%d",\ + fprintf(fichtm,"

      Parameter line 2

      • Tolerance for the convergence of the likelihood: ftol=%g \n
      • Interval for the elementary matrix (in month): stepm=%d",\ ftol, stepm); fprintf(fichtm,"\n
      • Number of fixed dummy covariates: ncovcol=%d ", ncovcol); ncurrv=1; @@ -11584,10 +11694,10 @@ Title=%s
        Datafile=%s Firstpass=%d La fprintf(fichtm,"\n
      • Number of fixed quantitative variables: nqv=%d ", nqv); ncurrv=i; for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i); - fprintf(fichtm,"\n
      • Number of time varying (wave varying) covariates: ntv=%d ", ntv); + fprintf(fichtm,"\n
      • Number of time varying (wave varying) dummy covariates: ntv=%d ", ntv); ncurrv=i; for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i); - fprintf(fichtm,"\n
      • Number of quantitative time varying covariates: nqtv=%d ", nqtv); + fprintf(fichtm,"\n
      • Number of time varying quantitative covariates: nqtv=%d ", nqtv); ncurrv=i; for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i); fprintf(fichtm,"\n
      • Weights column \n
        Number of alive states: nlstate=%d
        Number of death states (not really implemented): ndeath=%d \n
      • Number of waves: maxwav=%d \n
      • Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n
      • Does the weight column be taken into account (1), or not (0): weight=%d
      \n", \ @@ -11619,10 +11729,10 @@ Interval (in months) between two waves: for(j=1;j<=NDIM;j++) ximort[i][j]=0.; /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ - cens=ivector(1,n); - ageexmed=vector(1,n); - agecens=vector(1,n); - dcwave=ivector(1,n); + cens=ivector(firstobs,lastobs); + ageexmed=vector(firstobs,lastobs); + agecens=vector(firstobs,lastobs); + dcwave=ivector(firstobs,lastobs); for (i=1; i<=imx; i++){ dcwave[i]=-1; @@ -11836,9 +11946,10 @@ Please run with mle=-1 to get a correct free_vector(lpop,1,AGESUP); free_vector(tpop,1,AGESUP); free_matrix(ximort,1,NDIM,1,NDIM); - free_ivector(cens,1,n); - free_vector(agecens,1,n); - free_ivector(dcwave,1,n); + free_ivector(dcwave,firstobs,lastobs); + free_vector(agecens,firstobs,lastobs); + free_vector(ageexmed,firstobs,lastobs); + free_ivector(cens,firstobs,lastobs); #ifdef GSL #endif } /* Endof if mle==-3 mortality only */ @@ -11872,7 +11983,7 @@ Please run with mle=-1 to get a correct printf("\n"); /*--------- results files --------------*/ - 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,"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"); @@ -12244,8 +12355,8 @@ Please run with mle=-1 to get a correct /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ - free_lvector(num,1,n); - free_vector(agedc,1,n); + free_lvector(num,firstobs,lastobs); + free_vector(agedc,firstobs,lastobs); /*free_matrix(covar,0,NCOVMAX,1,n);*/ /*free_matrix(covar,1,NCOVMAX,1,n);*/ fclose(ficparo); @@ -12502,13 +12613,13 @@ Please run with mle=-1 to get a correct if(vpopbased==1) fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); else - fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); + fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n"); fprintf(ficrest,"# Age popbased mobilav e.. (std) "); for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ - printf("Computing age specific period (stable) prevalences in each health state \n"); - fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); + printf("Computing age specific forward period (stable) prevalences in each health state \n"); + fprintf(ficlog,"Computing age specific forward period (stable) prevalences in each health state \n"); for(age=bage; age <=fage ;age++){ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k, nres); /*ZZ Is it the correct prevalim */ if (vpopbased==1) { @@ -12555,16 +12666,16 @@ Please run with mle=-1 to get a correct printf("done State-specific expectancies\n");fflush(stdout); fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); - /* variance-covariance of period prevalence*/ + /* variance-covariance of forward period prevalence*/ varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); - free_vector(weight,1,n); + free_vector(weight,firstobs,lastobs); free_imatrix(Tvard,1,NCOVMAX,1,2); - free_imatrix(s,1,maxwav+1,1,n); - free_matrix(anint,1,maxwav,1,n); - free_matrix(mint,1,maxwav,1,n); - free_ivector(cod,1,n); + free_imatrix(s,1,maxwav+1,firstobs,lastobs); + free_matrix(anint,1,maxwav,firstobs,lastobs); + free_matrix(mint,1,maxwav,firstobs,lastobs); + free_ivector(cod,firstobs,lastobs); free_ivector(tab,1,NCOVMAX); fclose(ficresstdeij); fclose(ficrescveij); @@ -12584,10 +12695,10 @@ Please run with mle=-1 to get a correct 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); - if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n); - if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); - if(nqv>=1)free_matrix(coqvar,1,nqv,1,n); - free_matrix(covar,0,NCOVMAX,1,n); + if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs); + if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs); + if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs); + free_matrix(covar,0,NCOVMAX,firstobs,lastobs); free_matrix(matcov,1,npar,1,npar); free_matrix(hess,1,npar,1,npar); /*free_vector(delti,1,npar);*/ @@ -12740,5 +12851,6 @@ end: printf("\nType q for exiting: "); fflush(stdout); scanf("%s",z); } + printf("End\n"); exit(0); }