--- imach/src/imach.c 2019/06/04 13:51:20 1.301
+++ imach/src/imach.c 2022/04/05 21:03:51 1.311
@@ -1,6 +1,52 @@
-/* $Id: imach.c,v 1.301 2019/06/04 13:51:20 brouard Exp $
+/* $Id: imach.c,v 1.311 2022/04/05 21:03:51 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.311 2022/04/05 21:03:51 brouard
+ Summary: Fixed quantitative covariates
+
+ Fixed covariates (dummy or quantitative)
+ with missing values have never been allowed but are ERRORS and
+ program quits. Standard deviations of fixed covariates were
+ wrongly computed. Mean and standard deviations of time varying
+ covariates are still not computed.
+
+ Revision 1.310 2022/03/17 08:45:53 brouard
+ Summary: 99r25
+
+ Improving detection of errors: result lines should be compatible with
+ the model.
+
+ Revision 1.309 2021/05/20 12:39:14 brouard
+ Summary: Version 0.99r24
+
+ Revision 1.308 2021/03/31 13:11:57 brouard
+ Summary: Version 0.99r23
+
+
+ * imach.c (Module): Still bugs in the result loop. Thank to Holly Benett
+
+ Revision 1.307 2021/03/08 18:11:32 brouard
+ Summary: 0.99r22 fixed bug on result:
+
+ Revision 1.306 2021/02/20 15:44:02 brouard
+ Summary: Version 0.99r21
+
+ * imach.c (Module): Fix bug on quitting after result lines!
+ (Module): Version 0.99r21
+
+ Revision 1.305 2021/02/20 15:28:30 brouard
+ * imach.c (Module): Fix bug on quitting after result lines!
+
+ Revision 1.304 2021/02/12 11:34:20 brouard
+ * imach.c (Module): The use of a Windows BOM (huge) file is now an error
+
+ Revision 1.303 2021/02/11 19:50:15 brouard
+ * (Module): imach.c Someone entered 'results:' instead of 'result:'. Now it is an error which is printed.
+
+ Revision 1.302 2020/02/22 21:00:05 brouard
+ * (Module): imach.c Update mle=-3 (for computing Life expectancy
+ and life table from the data without any state)
+
Revision 1.301 2019/06/04 13:51:20 brouard
Summary: Error in 'r'parameter file backcast yearsbproj instead of yearsfproj
@@ -1122,12 +1168,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.301 2019/06/04 13:51:20 brouard Exp $ */
+/* $Id: imach.c,v 1.311 2022/04/05 21:03:51 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-char copyright[]="May 2019,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020";
-char fullversion[]="$Revision: 1.301 $ $Date: 2019/06/04 13:51:20 $";
+char copyright[]="March 2021,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021, INED 2000-2021";
+char fullversion[]="$Revision: 1.311 $ $Date: 2022/04/05 21:03:51 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1152,7 +1198,7 @@ int nqtveff=0; /**< ntqveff number of ef
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 npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */
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 */
@@ -1543,7 +1589,7 @@ char *cutl(char *blocc, char *alocc, cha
{
/* cuts string in into blocc and alocc where blocc ends before FIRST occurence of char 'occ'
and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2')
- gives blocc="abcdef" and alocc="ghi2j".
+ gives alocc="abcdef" and blocc="ghi2j".
If occ is not found blocc is null and alocc is equal to in. Returns blocc
*/
char *s, *t;
@@ -2415,13 +2461,13 @@ void powell(double p[], double **xi, int
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */
/* New value of last point Pn is not computed, P(n-1) */
for(j=1;j<=n;j++) {
- if(flatdir[j] >0){
- printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- }
- /* printf("\n"); */
- /* fprintf(ficlog,"\n"); */
- }
+ if(flatdir[j] >0){
+ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ }
+ /* printf("\n"); */
+ /* fprintf(ficlog,"\n"); */
+ }
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
@@ -3030,7 +3076,7 @@ double **pmij(double **ps, double *cov,
/* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */
double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij )
{
- /* Computes the backward probability at age agefin and covariate combination ij. In fact cov is already filled and x too.
+ /* Computes the backward probability at age agefin, cov[2], and covariate combination 'ij'. In fact cov is already filled and x too.
* Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij.
*/
int i, ii, j,k;
@@ -3810,7 +3856,7 @@ double funcone( double *x)
/* Fixed */
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
/* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
- for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */
+ for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
/* cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i]; */
/* cov[2+6]=covar[Tvar[6]][i]; */
@@ -4665,10 +4711,13 @@ 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 */
+ for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean on known values only */
+ if(!isnan(covar[ncovcol+z1][iind])){
+ 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]; *//*error*/
+ stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][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); */
@@ -4731,16 +4780,19 @@ Title=%s
Datafile=%s Firstpass=%d La
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]);
+ fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.3g (weighted) individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]);
+ fprintf(ficlog,", mean=%.3g",meanq[z1]/idq[z1]," stdeviation=%.3g\n",stdq[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])); + /* mu = \frac{w x}{\sum w} + var = \frac{\sum w (x-mu)^2}{\sum w} = \frac{w x^2}{\sum w} - mu^2 + */ + printf(" fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); + fprintf(ficlog," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); + fprintf(ficresphtmfr," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)
\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); } /* for (z1=1; z1<= nqtveff; z1++) { */ /* for(m=1;m<=lastpass;m++){ */ @@ -5243,33 +5295,33 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ mi=0; /* First valid wave */ mli=0; /* Last valid wave */ - m=firstpass; - while(s[m][i] <= nlstate){ /* a live state */ + m=firstpass; /* Loop on waves */ + while(s[m][i] <= nlstate){ /* a live state or unknown state */ if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */ mli=m-1;/* mw[++mi][i]=m-1; */ }else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ - mw[++mi][i]=m; + mw[++mi][i]=m; /* Valid wave: incrementing mi and updating mi; mw[mi] is the wave number of mi_th valid transition */ mli=m; } /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */ if(m < lastpass){ /* m < lastpass, standard case */ m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */ } - else{ /* m >= lastpass, eventual special issue with warning */ + else{ /* m = lastpass, eventual special issue with warning */ #ifdef UNKNOWNSTATUSNOTCONTRIBUTING break; #else - if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ + if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ /* case -2 (vital status unknown is warned later */ if(firsthree == 0){ - printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\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, s[m][i], nlstate+ndeath); + printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\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, s[m][i], nlstate+ndeath); 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 as 1-p%d%d .\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, s[m][i], nlstate+ndeath); - mw[++mi][i]=m; + fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\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, s[m][i], nlstate+ndeath); + mw[++mi][i]=m; /* Valid transition with unknown status */ mli=m; } if(s[m][i]==-2){ /* Vital status is really unknown */ nbwarn++; - if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ + if((int)anint[m][i] == 9999){ /* Has the vital status really been verified?not a transition */ 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.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); } @@ -5294,34 +5346,35 @@ void concatwav(int wav[], int **dh, int #ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE else if ((int) andc[i] != 9999) { /* Date of death is known */ if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ - if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */ + if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* month of death occured before last wave month and status should have been death instead of -1 */ nbwarn++; if(firstfiv==0){ - printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d, interviewed on %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); firstfiv=1; }else{ - fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d, interviewed on %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); } - }else{ /* Death occured afer last wave potential bias */ + s[m][i]=nlstate+1; /* Fixing the status as death. Be careful if multiple death states */ + }else{ /* Month of Death occured afer last wave month, potential bias */ nberr++; if(firstwo==0){ - printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\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 ); + printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d with status %d. Potential bias if other individuals are still alive on this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictitious wave at the date of last vital status scan, with a dead status. See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); firstwo=1; } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d with status %d. Potential bias if other individuals are still alive on this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictitious wave at the date of last vital status scan, with a dead status. See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); } }else{ /* if date of interview is unknown */ /* 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 ); + 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 with status %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], s[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 ); + 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 with status %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], s[m][i], i,m ); } } /* end if date of death is known */ #endif - wav[i]=mi; /* mi should be the last effective wave (or mli) */ - /* wav[i]=mw[mi][i]; */ + wav[i]=mi; /* mi should be the last effective wave (or mli), */ + /* wav[i]=mw[mi][i]; */ if(mi==0){ nbwarn++; if(first==0){ @@ -5456,6 +5509,8 @@ void concatwav(int wav[], int **dh, int if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ switch(Fixed[k]) { case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ + modmaxcovj=0; + modmincovj=0; for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ ij=(int)(covar[Tvar[k]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i @@ -5469,8 +5524,10 @@ void concatwav(int wav[], int **dh, int else if (ij < modmincovj) modmincovj=ij; 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); + printf("ERROR, IMaCh doesn't treat covariate with missing values V%d=-1, individual %d will be skipped.\n",Tvar[k],i); + fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=-1, individual %d will be skipped.\n",Tvar[k],i); + fflush(ficlog); + exit(1); } if ((ij < -1) || (ij > NCOVMAX)){ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); @@ -5545,6 +5602,16 @@ void concatwav(int wav[], int **dh, int break; } /* end switch */ } /* end dummy test */ + if(Dummy[k]==1 && Typevar[k] !=1){ /* Dummy covariate and not age product */ + for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ + if(isnan(covar[Tvar[k]][i])){ + printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); + fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); + fflush(ficlog); + exit(1); + } + } + } } /* 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; @@ -8505,7 +8572,7 @@ void prevforecast(char fileres[], double */ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; double agec; /* generic age */ - double agelim, ppij, ppi, yp,yp1,yp2,jintmean,mintmean,aintmean; + double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/ double *popeffectif,*popcount; double ***p3mat; /* double ***mobaverage; */ @@ -9061,7 +9128,7 @@ void prwizard(int ncovmodel, int nlstate /******************* Gompertz Likelihood ******************************/ double gompertz(double x[]) { - double A,B,L=0.0,sump=0.,num=0.; + double A=0.0,B=0.,L=0.0,sump=0.,num=0.; int i,n=0; /* n is the size of the sample */ for (i=1;i<=imx ; i++) { @@ -9069,28 +9136,34 @@ double gompertz(double x[]) /* sump=sump+1;*/ num=num+1; } - - + L=0.0; + /* agegomp=AGEGOMP; */ /* for (i=0; i<=imx; i++) if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ - for (i=1;i<=imx ; i++) - { - if (cens[i] == 1 && wav[i]>1) - A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); - - if (cens[i] == 0 && wav[i]>1) + for (i=1;i<=imx ; i++) { + /* mu(a)=mu(agecomp)*exp(teta*(age-agegomp)) + mu(a)=x[1]*exp(x[2]*(age-agegomp)); x[1] and x[2] are per year. + * L= Product mu(agedeces)exp(-\int_ageexam^agedc mu(u) du ) for a death between agedc (in month) + * and agedc +1 month, cens[i]=0: log(x[1]/YEARM) + * + + * exp(-\int_ageexam^agecens mu(u) du ) when censored, cens[i]=1 + */ + if (wav[i] > 1 || agedc[i] < AGESUP) { + if (cens[i] == 1){ + A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); + } else if (cens[i] == 0){ A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) - +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); - + +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); + } else + printf("Gompertz cens[%d] neither 1 nor 0\n",i); /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ - if (wav[i] > 1 ) { /* ??? */ - L=L+A*weight[i]; + L=L+A*weight[i]; /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ - } - } + } + } - /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ + /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ return -2*L*num/sump; } @@ -9099,7 +9172,7 @@ double gompertz(double x[]) /******************* Gompertz_f Likelihood ******************************/ double gompertz_f(const gsl_vector *v, void *params) { - double A,B,LL=0.0,sump=0.,num=0.; + double A=0.,B=0.,LL=0.0,sump=0.,num=0.; double *x= (double *) v->data; int i,n=0; /* n is the size of the sample */ @@ -9192,6 +9265,7 @@ int readdata(char datafile[], int firsto int i=0, j=0, n=0, iv=0, v; int lstra; int linei, month, year,iout; + int noffset=0; /* This is the offset if BOM data file */ char line[MAXLINE], linetmp[MAXLINE]; char stra[MAXLINE], strb[MAXLINE]; char *stratrunc; @@ -9225,8 +9299,53 @@ int readdata(char datafile[], int firsto fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1; } - i=1; + /* Is it a BOM UTF-8 Windows file? */ + /* First data line */ linei=0; + while(fgets(line, MAXLINE, fic)) { + noffset=0; + if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */ + { + noffset=noffset+3; + printf("# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile);fflush(stdout); + fprintf(ficlog,"# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile); + fflush(ficlog); return 1; + } + /* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/ + else if( line[0] == (char)0xFF && line[1] == (char)0xFE) + { + noffset=noffset+2; + printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); + fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); + fflush(ficlog); return 1; + } + else if( line[0] == 0 && line[1] == 0) + { + if( line[2] == (char)0xFE && line[3] == (char)0xFF){ + noffset=noffset+4; + printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); + fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); + fflush(ficlog); return 1; + } + } else{ + ;/*printf(" Not a BOM file\n");*/ + } + /* If line starts with a # it is a comment */ + if (line[noffset] == '#') { + linei=linei+1; + break; + }else{ + break; + } + } + fclose(fic); + if((fic=fopen(datafile,"r"))==NULL) { + printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout); + fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1; + } + /* Not a Bom file */ + + i=1; while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { linei=linei+1; for(j=strlen(line); j>=0;j--){ /* Untabifies line */ @@ -9347,7 +9466,11 @@ int readdata(char datafile[], int firsto return 1; } anint[j][i]= (double) year; - mint[j][i]= (double)month; + mint[j][i]= (double)month; + /* if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ */ + /* printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */ + /* fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */ + /* } */ strcpy(line,stra); } /* End loop on waves */ @@ -9386,7 +9509,14 @@ int readdata(char datafile[], int firsto } annais[i]=(double)(year); - moisnais[i]=(double)(month); + moisnais[i]=(double)(month); + for (j=1;j<=maxwav;j++){ + if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ + printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j,(int)moisnais[i],(int)annais[i]); + fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j, (int)moisnais[i],(int)annais[i]); + } + } + strcpy(line,stra); /* Sample weight */ @@ -9406,6 +9536,8 @@ int readdata(char datafile[], int firsto cutv(stra, strb, line, ' '); if(strb[0]=='.') { /* Missing value */ lval=-1; + coqvar[iv][i]=NAN; + covar[ncovcol+iv][i]=NAN; /* including qvar in standard covar for performance reasons */ }else{ errno=0; /* what_kind_of_number(strb); */ @@ -9509,7 +9641,6 @@ int decoderesult ( char resultline[], in char stra[80], strb[80], strc[80], strd[80],stre[80]; removefirstspace(&resultline); - printf("decoderesult:%s\n",resultline); if (strstr(resultline,"v") !=0){ printf("Error. 'v' must be in upper case 'V' result: %s ",resultline); @@ -9524,15 +9655,14 @@ int decoderesult ( char resultline[], in TKresult[nres]=0; /* Combination for the nresult and the model */ return (0); } - if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ - printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); - fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); + printf("ERROR: the number of variables in the resultline, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); + fprintf(ficlog,"ERROR: the number of variables in the resultline, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); } for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */ if(nbocc(resultsav,'=') >1){ cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' - resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */ + resultsav= V4=1 V5=25.1 V3=0 stra= V5=25.1 V3=0 strb= V4=1 */ cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */ }else cutl(strc,strd,resultsav,'='); @@ -9557,7 +9687,9 @@ int decoderesult ( char resultline[], in } } if(match == 0){ - printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + printf("Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); + fprintf(ficlog,"Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); + return 1; } } } @@ -9574,8 +9706,12 @@ int decoderesult ( char resultline[], in } if(match == 0){ printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + fprintf(ficlog,"Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + return 1; }else if(match > 1){ printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); + fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); + return 1; } } @@ -11085,7 +11221,8 @@ int main(int argc, char *argv[]) noffset=noffset+3; printf("# File is an UTF8 Bom.\n"); // 0xBF } - else if( line[0] == (char)0xFE && line[1] == (char)0xFF) +/* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/ + else if( line[0] == (char)0xFF && line[1] == (char)0xFE) { noffset=noffset+2; printf("# File is an UTF16BE BOM file\n"); @@ -11173,8 +11310,8 @@ int main(int argc, char *argv[]) } if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ if (num_filled != 1){ - printf("ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line); - fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line); + printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); + fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); model[0]='\0'; goto end; } @@ -11860,8 +11997,8 @@ Interval (in months) between two waves: ximort[i][j]=(i == j ? 1.0 : 0.0); } - /*p[1]=0.0268; p[NDIM]=0.083;*/ - /*printf("%lf %lf", p[1], p[2]);*/ + p[1]=0.0268; p[NDIM]=0.083; + /* printf("%lf %lf", p[1], p[2]); */ #ifdef GSL @@ -11987,9 +12124,9 @@ Interval (in months) between two waves: printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); } - lsurv=vector(1,AGESUP); - lpop=vector(1,AGESUP); - tpop=vector(1,AGESUP); + lsurv=vector(agegomp,AGESUP); + lpop=vector(agegomp,AGESUP); + tpop=vector(agegomp,AGESUP); lsurv[agegomp]=100000; for (k=agegomp;k<=AGESUP;k++) { @@ -12036,9 +12173,9 @@ Please run with mle=-1 to get a correct stepm, weightopt,\ model,imx,p,matcov,agemortsup); - free_vector(lsurv,1,AGESUP); - free_vector(lpop,1,AGESUP); - free_vector(tpop,1,AGESUP); + free_vector(lsurv,agegomp,AGESUP); + free_vector(lpop,agegomp,AGESUP); + free_vector(tpop,agegomp,AGESUP); free_matrix(ximort,1,NDIM,1,NDIM); free_ivector(dcwave,firstobs,lastobs); free_vector(agecens,firstobs,lastobs); @@ -12334,11 +12471,13 @@ Please run with mle=-1 to get a correct } /* Results */ + endishere=0; nresult=0; + parameterline=0; do{ if(!fgets(line, MAXLINE, ficpar)){ endishere=1; - parameterline=14; + parameterline=15; }else if (line[0] == '#') { /* If line starts with a # it is a comment */ numlinepar++; @@ -12351,12 +12490,13 @@ Please run with mle=-1 to get a correct parameterline=11; else if(sscanf(line,"prevbackcast=%[^\n]\n",modeltemp)) parameterline=12; - else if(sscanf(line,"result:%[^\n]\n",modeltemp)) + else if(sscanf(line,"result:%[^\n]\n",modeltemp)){ parameterline=13; + } else{ parameterline=14; } - switch (parameterline){ + switch (parameterline){ /* =0 only if only comments */ case 11: if((num_filled=sscanf(line,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj)) !=EOF && (num_filled == 8)){ fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); @@ -12369,9 +12509,9 @@ Please run with mle=-1 to get a correct prvforecast = 1; } else if((num_filled=sscanf(line,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",&prevfcast,&yrfproj,&mobilavproj)) !=EOF){/* && (num_filled == 3))*/ - printf("prevforecast=%d yearsfproj=%lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); - fprintf(ficlog,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); - fprintf(ficres,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); + printf("prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); + fprintf(ficlog,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); + fprintf(ficres,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); prvforecast = 2; } else { @@ -12392,9 +12532,9 @@ Please run with mle=-1 to get a correct prvbackcast = 1; } else if((num_filled=sscanf(line,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",&prevbcast,&yrbproj,&mobilavproj)) ==3){/* && (num_filled == 3))*/ - printf("prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); - fprintf(ficlog,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); - fprintf(ficres,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); + printf("prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); + fprintf(ficlog,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); + fprintf(ficres,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); prvbackcast = 2; } else { @@ -12404,38 +12544,37 @@ Please run with mle=-1 to get a correct } break; case 13: - if((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){ - if (num_filled == 0){ - resultline[0]='\0'; - printf("Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line); - fprintf(ficlog,"Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line); - break; - } else if (num_filled != 1){ - printf("ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line); - fprintf(ficlog,"ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line); - } - nresult++; /* Sum of resultlines */ - printf("Result %d: result=%s\n",nresult, resultline); - if(nresult > MAXRESULTLINES){ - printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult); - fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult); - goto end; - } - decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */ - fprintf(ficparo,"result: %s\n",resultline); - fprintf(ficres,"result: %s\n",resultline); - fprintf(ficlog,"result: %s\n",resultline); - break; - case 14: - if(ncovmodel >2 && nresult==0 ){ - printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); - goto end; - } - break; - default: - nresult=1; - decoderesult(".",nresult ); /* No covariate */ + num_filled=sscanf(line,"result:%[^\n]\n",resultline); + nresult++; /* Sum of resultlines */ + printf("Result %d: result:%s\n",nresult, resultline); + if(nresult > MAXRESULTLINES){ + printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINES,nresult,rfileres); + fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINES,nresult,rfileres); + goto end; } + if(!decoderesult(resultline, nresult)){ /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */ + fprintf(ficparo,"result: %s\n",resultline); + fprintf(ficres,"result: %s\n",resultline); + fprintf(ficlog,"result: %s\n",resultline); + } else + goto end; + break; + case 14: + printf("Error: Unknown command '%s'\n",line); + fprintf(ficlog,"Error: Unknown command '%s'\n",line); + if(ncovmodel >=2 && nresult==0 ){ + printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); + fprintf(ficlog,"ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); + } + /* goto end; */ + break; + case 15: + printf("End of resultlines.\n"); + fprintf(ficlog,"End of resultlines.\n"); + break; + default: /* parameterline =0 */ + nresult=1; + decoderesult(".",nresult ); /* No covariate */ } /* End switch parameterline */ }while(endishere==0); /* End do */