--- imach/src/imach.c 2016/08/22 06:55:53 1.230 +++ imach/src/imach.c 2016/08/22 07:17:15 1.231 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.230 2016/08/22 06:55:53 brouard Exp $ +/* $Id: imach.c,v 1.231 2016/08/22 07:17:15 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.231 2016/08/22 07:17:15 brouard + Summary: not working + Revision 1.230 2016/08/22 06:55:53 brouard Summary: Not working @@ -890,12 +893,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.230 2016/08/22 06:55:53 brouard Exp $ */ +/* $Id: imach.c,v 1.231 2016/08/22 07:17:15 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.230 $ $Date: 2016/08/22 06:55:53 $"; +char fullversion[]="$Revision: 1.231 $ $Date: 2016/08/22 07:17:15 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1064,6 +1067,15 @@ double ***cotvar; /* Time varying covari double ***cotqvar; /* Time varying quantitative covariate itqv */ double idx; int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ +int *TvarFD; /**< TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +int *TvarFDind; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +int *TvarFQ; /* TvarFQ[1]=V2 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ +int *TvarFQind; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ +int *TvarVD; /* TvarVD[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ +int *TvarVDind; /* TvarVDind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ +int *TvarVQ; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ +int *TvarVQind; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ + int *Tvarsel; /**< Selected covariates for output */ double *Tvalsel; /**< Selected modality value of covariate for output */ int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ @@ -1086,6 +1098,33 @@ int *Tposprod; /**< Gives the k1 product int cptcovprod, *Tvaraff, *invalidvarcomb; double *lsurv, *lpop, *tpop; +#define FD 1; /* Fixed dummy covariate */ +#define FQ 2; /* Fixed quantitative covariate */ +#define FP 3; /* Fixed product covariate */ +#define FPDD 7; /* Fixed product dummy*dummy covariate */ +#define FPDQ 8; /* Fixed product dummy*quantitative covariate */ +#define FPQQ 9; /* Fixed product quantitative*quantitative covariate */ +#define VD 10; /* Varying dummy covariate */ +#define VQ 11; /* Varying quantitative covariate */ +#define VP 12; /* Varying product covariate */ +#define VPDD 13; /* Varying product dummy*dummy covariate */ +#define VPDQ 14; /* Varying product dummy*quantitative covariate */ +#define VPQQ 15; /* Varying product quantitative*quantitative covariate */ +#define APFD 16; /* Age product * fixed dummy covariate */ +#define APFQ 17; /* Age product * fixed quantitative covariate */ +#define APVD 18; /* Age product * varying dummy covariate */ +#define APVQ 19; /* Age product * varying quantitative covariate */ + +#define FTYPE 1; /* Fixed covariate */ +#define VTYPE 2; /* Varying covariate (loop in wave) */ +#define ATYPE 2; /* Age product covariate (loop in dh within wave)*/ + +struct kmodel{ + int maintype; /* main type */ + int subtype; /* subtype */ +}; +struct kmodel modell[NCOVMAX]; + double ftol=FTOL; /**< Tolerance for computing Max Likelihood */ double ftolhess; /**< Tolerance for computing hessian */ @@ -3006,94 +3045,94 @@ double func( double *x) meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] */ for(mi=1; mi<= wav[i]-1; mi++){ - for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */ - /* cov[ioffset+itv]=cotvar[mw[mi][i]][Tvar[itv]][i]; /\* Not sure, Tvar V4+V3+V5 Tvaraff ? *\/ */ - cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; - } - for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */ - if(cotqvar[mw[mi][i]][iqtv][i] == -1){ - printf("i=%d, mi=%d, iqtv=%d, cotqvar[mw[mi][i]][iqtv][i]=%f",i,mi,iqtv,cotqvar[mw[mi][i]][iqtv][i]); - } - cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; - /* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i]; */ - } - /* ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; */ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } - for(d=0; d 1 the results are less biased than in previous versions. - */ - s1=s[mw[mi][i]][i]; - s2=s[mw[mi+1][i]][i]; - bbh=(double)bh[mi][i]/(double)stepm; - /* bias bh is positive if real duration - * is higher than the multiple of stepm and negative otherwise. - */ - /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ - if( s2 > nlstate){ - /* i.e. if s2 is a death state and if the date of death is known - then the contribution to the likelihood is the probability to - die between last step unit time and current step unit time, - which is also equal to probability to die before dh - minus probability to die before dh-stepm . - In version up to 0.92 likelihood was computed - as if date of death was unknown. Death was treated as any other - health state: the date of the interview describes the actual state - and not the date of a change in health state. The former idea was - to consider that at each interview the state was recorded - (healthy, disable or death) and IMaCh was corrected; but when we - introduced the exact date of death then we should have modified - the contribution of an exact death to the likelihood. This new - contribution is smaller and very dependent of the step unit - stepm. It is no more the probability to die between last interview - and month of death but the probability to survive from last - interview up to one month before death multiplied by the - probability to die within a month. Thanks to Chris - Jackson for correcting this bug. Former versions increased - mortality artificially. The bad side is that we add another loop - which slows down the processing. The difference can be up to 10% - lower mortality. - */ - /* If, at the beginning of the maximization mostly, the - cumulative probability or probability to be dead is - constant (ie = 1) over time d, the difference is equal to - 0. out[s1][3] = savm[s1][3]: probability, being at state - s1 at precedent wave, to be dead a month before current - wave is equal to probability, being at state s1 at - precedent wave, to be dead at mont of the current - wave. Then the observed probability (that this person died) - is null according to current estimated parameter. In fact, - it should be very low but not zero otherwise the log go to - infinity. - */ + /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ + /* But now since version 0.9 we anticipate for bias at large stepm. + * If stepm is larger than one month (smallest stepm) and if the exact delay + * (in months) between two waves is not a multiple of stepm, we rounded to + * the nearest (and in case of equal distance, to the lowest) interval but now + * we keep into memory the bias bh[mi][i] and also the previous matrix product + * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the + * probability in order to take into account the bias as a fraction of the way + * from savm to out if bh is negative or even beyond if bh is positive. bh varies + * -stepm/2 to stepm/2 . + * For stepm=1 the results are the same as for previous versions of Imach. + * For stepm > 1 the results are less biased than in previous versions. + */ + s1=s[mw[mi][i]][i]; + s2=s[mw[mi+1][i]][i]; + bbh=(double)bh[mi][i]/(double)stepm; + /* bias bh is positive if real duration + * is higher than the multiple of stepm and negative otherwise. + */ + /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ + if( s2 > nlstate){ + /* i.e. if s2 is a death state and if the date of death is known + then the contribution to the likelihood is the probability to + die between last step unit time and current step unit time, + which is also equal to probability to die before dh + minus probability to die before dh-stepm . + In version up to 0.92 likelihood was computed + as if date of death was unknown. Death was treated as any other + health state: the date of the interview describes the actual state + and not the date of a change in health state. The former idea was + to consider that at each interview the state was recorded + (healthy, disable or death) and IMaCh was corrected; but when we + introduced the exact date of death then we should have modified + the contribution of an exact death to the likelihood. This new + contribution is smaller and very dependent of the step unit + stepm. It is no more the probability to die between last interview + and month of death but the probability to survive from last + interview up to one month before death multiplied by the + probability to die within a month. Thanks to Chris + Jackson for correcting this bug. Former versions increased + mortality artificially. The bad side is that we add another loop + which slows down the processing. The difference can be up to 10% + lower mortality. + */ + /* If, at the beginning of the maximization mostly, the + cumulative probability or probability to be dead is + constant (ie = 1) over time d, the difference is equal to + 0. out[s1][3] = savm[s1][3]: probability, being at state + s1 at precedent wave, to be dead a month before current + wave is equal to probability, being at state s1 at + precedent wave, to be dead at mont of the current + wave. Then the observed probability (that this person died) + is null according to current estimated parameter. In fact, + it should be very low but not zero otherwise the log go to + infinity. + */ /* #ifdef INFINITYORIGINAL */ /* lli=log(out[s1][s2] - savm[s1][s2]); */ /* #else */ @@ -3314,53 +3353,56 @@ double funcone( double *x) for (i=1,ipmx=0, sw=0.; i<=imx; i++){ ioffset=2+nagesqr+cptcovage; /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ - for (k=1; k<=ncoveff+nqfveff;k++){ /* Simple and product fixed Dummy covariates without age* products */ + for (k=1; k<=ncoveff;k++){ /* Simple and product fixed Dummy covariates without age* products */ cov[++ioffset]=covar[TvarFD[k]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ } - for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitative fixed covariates */ - cov[++ioffset]=coqvar[Tvar[iqv]][i]; /* Only V1 k=9 */ + for (k=1; k<=nqfveff;k++){ /* Simple and product fixed Quantitative covariates without age* products */ + cov[++ioffset]=coqvar[TvarFQ[k]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V2 and V1*V2 is fixed (k=6 and 7?)*/ } + /* for(iqv=1; iqv <= nqfveff; iqv++){ /\* Quantitative fixed covariates *\/ */ + /* cov[++ioffset]=coqvar[Tvar[iqv]][i]; /\* Only V2 k=6 and V1*V2 7 *\/ */ + /* } */ for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ - for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */ - /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ - /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ - k=ioffset-2-nagesqr-cptcovage+itv; /* position in simple model */ - cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; - /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ + for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates (single??)*/ + /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ + /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ + k=ioffset-2-nagesqr-cptcovage+itv; /* position in simple model */ + cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; + /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ } for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */ - iv=TmodelInvQind[iqtv]; /* Counting the # varying covariate from 1 to ntveff */ - /* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); */ - cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; + iv=TmodelInvQind[iqtv]; /* Counting the # varying covariate from 1 to ntveff */ + /* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); */ + cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; } for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ for(d=0; dDatafile=%s Firstpass=%d La scanf("%d", i);*/ for (i=-5; i<=nlstate+ndeath; i++) for (jk=-5; jk<=nlstate+ndeath; jk++) - for(m=iagemin; m <= iagemax+3; m++) - freq[i][jk][m]=0; - + for(m=iagemin; m <= iagemax+3; m++) + freq[i][jk][m]=0; + for (i=1; i<=nlstate; i++) { for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0; + prop[i][m]=0; posprop[i]=0; pospropt[i]=0; } @@ -4011,87 +4053,87 @@ Title=%s
Datafile=%s Firstpass=%d La /* meanqt[m][z1]=0.; */ /* } */ /* } */ - + dateintsum=0; k2cpt=0; /* For that combination of covariate j1, we count and print the frequencies in one pass */ for (iind=1; iind<=imx; iind++) { /* For each individual iind */ bool=1; if(anyvaryingduminmodel==0){ /* If All fixed covariates */ - if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + 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++) { - /* if(Tvaraff[z1] ==-20){ */ - /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ - /* }else if(Tvaraff[z1] ==-10){ */ - /* /\* sumnew+=coqvar[z1][iind]; *\/ */ - /* }else */ - if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ - /* Tests if this individual iind responded to j1 (V4=1 V3=0) */ - bool=0; - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", - bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), - j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ - /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ - } /* Onlyf fixed */ - } /* end z1 */ - } /* cptcovn > 0 */ + for (z1=1; z1<=cptcoveff; z1++) { + /* if(Tvaraff[z1] ==-20){ */ + /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ + /* }else if(Tvaraff[z1] ==-10){ */ + /* /\* sumnew+=coqvar[z1][iind]; *\/ */ + /* }else */ + if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ + /* Tests if this individual iind responded to j1 (V4=1 V3=0) */ + bool=0; + /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", + bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), + j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ + /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ + } /* Onlyf fixed */ + } /* end z1 */ + } /* cptcovn > 0 */ } /* end any */ if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */ - /* for(m=firstpass; m<=lastpass; m++){ */ - for(mi=1; mi=firstpass && m <=lastpass){ - k2=anint[m][iind]+(mint[m][iind]/12.); - /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ - if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ - if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ - if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ - prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ - if (m1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { - dateintsum=dateintsum+k2; - k2cpt++; - /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ - } - } /* end bool 2 */ - } /* end m */ + /* for(m=firstpass; m<=lastpass; m++){ */ + for(mi=1; mi=firstpass && m <=lastpass){ + k2=anint[m][iind]+(mint[m][iind]/12.); + /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ + if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ + if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ + if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ + prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ + if (m1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { + dateintsum=dateintsum+k2; + k2cpt++; + /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ + } + } /* end bool 2 */ + } /* end m */ } /* end bool */ } /* end iind = 1 to imx */ /* prop[s][age] is feeded for any initial and valid live state as well as freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ - - + + /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ pstamp(ficresp); /* if (ncoveff>0) { */ @@ -4100,9 +4142,9 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtm, "\n

********** Variable "); fprintf(ficresphtmfr, "\n

********** Variable "); for (z1=1; z1<=cptcoveff; z1++){ - fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); } fprintf(ficresp, "**********\n#"); fprintf(ficresphtm, "**********

\n"); @@ -4118,126 +4160,126 @@ Title=%s
Datafile=%s Firstpass=%d La } fprintf(ficresp, "\n"); fprintf(ficresphtm, "\n"); - + /* Header of frequency table by age */ fprintf(ficresphtmfr,""); fprintf(ficresphtmfr," "); for(jk=-1; jk <=nlstate+ndeath; jk++){ for(m=-1; m <=nlstate+ndeath; m++){ - if(jk!=0 && m!=0) - fprintf(ficresphtmfr," ",jk,m); + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",jk,m); } } fprintf(ficresphtmfr, "\n"); - + /* For each age */ for(iage=iagemin; iage <= iagemax+3; iage++){ fprintf(ficresphtm,""); if(iage==iagemax+1){ - fprintf(ficlog,"1"); - fprintf(ficresphtmfr," "); + fprintf(ficlog,"1"); + fprintf(ficresphtmfr," "); }else if(iage==iagemax+2){ - fprintf(ficlog,"0"); - fprintf(ficresphtmfr," "); + fprintf(ficlog,"0"); + fprintf(ficresphtmfr," "); }else if(iage==iagemax+3){ - fprintf(ficlog,"Total"); - fprintf(ficresphtmfr," "); + fprintf(ficlog,"Total"); + fprintf(ficresphtmfr," "); }else{ - if(first==1){ - first=0; - printf("See log file for details...\n"); - } - fprintf(ficresphtmfr," ",iage); - fprintf(ficlog,"Age %d", iage); + if(first==1){ + first=0; + printf("See log file for details...\n"); + } + fprintf(ficresphtmfr," ",iage); + fprintf(ficlog,"Age %d", iage); } for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][iage]; + for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) + pp[jk] += freq[jk][m][iage]; } for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][iage]; - if(pp[jk]>=1.e-10){ - if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); - } - fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); - }else{ - if(first==1) - printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - } + for(m=-1, pos=0; m <=0 ; m++) + pos += freq[jk][m][iage]; + if(pp[jk]>=1.e-10){ + if(first==1){ + printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + } + fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + }else{ + if(first==1) + printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); + fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); + } } - + for(jk=1; jk <=nlstate ; jk++){ - /* posprop[jk]=0; */ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ - pp[jk] += freq[jk][m][iage]; + /* posprop[jk]=0; */ + for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ + pp[jk] += freq[jk][m][iage]; } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ - + for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ - pos += pp[jk]; /* pos is the total number of transitions until this age */ - posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ - pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pos += pp[jk]; /* pos is the total number of transitions until this age */ + posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state + from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state + from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ } for(jk=1; jk <=nlstate ; jk++){ - if(pos>=1.e-5){ - if(first==1) - printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - }else{ - if(first==1) - printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - } - if( iage <= iagemax){ - if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - /*probs[iage][jk][j1]= pp[jk]/pos;*/ - /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ - } - else{ - fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage, prop[jk][iage],pospropta); - } - } - pospropt[jk] +=posprop[jk]; + if(pos>=1.e-5){ + if(first==1) + printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + }else{ + if(first==1) + printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); + fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); + } + if( iage <= iagemax){ + if(pos>=1.e-5){ + fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); + fprintf(ficresphtm,"",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); + /*probs[iage][jk][j1]= pp[jk]/pos;*/ + /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ + } + else{ + fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); + fprintf(ficresphtm,"",iage, prop[jk][iage],pospropta); + } + } + pospropt[jk] +=posprop[jk]; } /* end loop jk */ /* pospropt=0.; */ for(jk=-1; jk <=nlstate+ndeath; jk++){ - for(m=-1; m <=nlstate+ndeath; m++){ - if(freq[jk][m][iage] !=0 ) { /* minimizing output */ - if(first==1){ - printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); - } - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); - } - if(jk!=0 && m!=0) - fprintf(ficresphtmfr," ",freq[jk][m][iage]); - } + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[jk][m][iage] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); + } + fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); + } + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",freq[jk][m][iage]); + } } /* end loop jk */ posproptt=0.; for(jk=1; jk <=nlstate; jk++){ - posproptt += pospropt[jk]; + posproptt += pospropt[jk]; } fprintf(ficresphtmfr,"\n "); if(iage <= iagemax){ - fprintf(ficresp,"\n"); - fprintf(ficresphtm,"\n"); + fprintf(ficresp,"\n"); + fprintf(ficresphtm,"\n"); } if(first==1) - printf("Others in log...\n"); + printf("Others in log...\n"); fprintf(ficlog,"\n"); } /* end loop age iage */ fprintf(ficresphtm,""); for(jk=1; jk <=nlstate ; jk++){ if(posproptt < 1.e-5){ - fprintf(ficresphtm,"",pospropt[jk],posproptt); + fprintf(ficresphtm,"",pospropt[jk],posproptt); }else{ - fprintf(ficresphtm,"",pospropt[jk]/posproptt,pospropt[jk],posproptt); + fprintf(ficresphtm,"",pospropt[jk]/posproptt,pospropt[jk],posproptt); } } fprintf(ficresphtm,"\n"); @@ -4255,7 +4297,7 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtmfr,"
Age%d%d%d%d
0
0
Unknown
Unknown
Total
Total
%d
%d%d%.5f%.0f%.0f%dNaNq%.0f%.0f%d%.5f%.0f%.0f%dNaNq%.0f%.0f%.0f%.0f
TotNanq%.0f%.0fNanq%.0f%.0f%.5f%.0f%.0f%.5f%.0f%.0f
\n"); } /* end selected combination of covariate j1 */ dateintmean=dateintsum/k2cpt; - + fclose(ficresp); fclose(ficresphtm); fclose(ficresphtmfr); @@ -4616,85 +4658,83 @@ 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 */ - 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 - * If product of Vn*Vm, still boolean *: - * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables - * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */ - /* Finds for covariate j, n=Tvar[j] of Vn . ij is the - modality of the nth covariate of individual i. */ - if (ij > modmaxcovj) - modmaxcovj=ij; - else if (ij < modmincovj) - modmincovj=ij; - if ((ij < -1) && (ij > NCOVMAX)){ - printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); - exit(1); - }else - Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/ - /* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */ - /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ - /* getting the maximum value of the modality of the covariate - (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and - female ies 1, then modmaxcovj=1. - */ - } /* end for loop on individuals i */ - printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); - fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); - cptcode=modmaxcovj; - /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ - /*for (i=0; i<=cptcode; i++) {*/ - for (j=modmincovj; j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */ - printf("Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); - fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); - if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */ - if( j != -1){ - ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th - covariate for which somebody answered excluding - undefined. Usually 2: 0 and 1. */ - } - ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th - covariate for which somebody answered including - undefined. Usually 3: -1, 0 and 1. */ - } - /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for - * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ - } /* Ndum[-1] number of undefined modalities */ - - /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ - /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. - If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; - modmincovj=3; modmaxcovj = 7; - There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; - which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; - defining two dummy variables: variables V1_1 and V1_2. - nbcode[Tvar[j]][ij]=k; - nbcode[Tvar[j]][1]=0; - nbcode[Tvar[j]][2]=1; - 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*/ - 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*/ - cptcode = ij; /* New max modality for covar j */ - } /* end of loop on modality i=-1 to 1 or more */ - break; + 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 + * If product of Vn*Vm, still boolean *: + * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables + * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */ + /* Finds for covariate j, n=Tvar[j] of Vn . ij is the + modality of the nth covariate of individual i. */ + if (ij > modmaxcovj) + modmaxcovj=ij; + else if (ij < modmincovj) + modmincovj=ij; + if ((ij < -1) && (ij > NCOVMAX)){ + printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); + exit(1); + }else + Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/ + /* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */ + /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ + /* getting the maximum value of the modality of the covariate + (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and + female ies 1, then modmaxcovj=1. + */ + } /* end for loop on individuals i */ + printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); + fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); + cptcode=modmaxcovj; + /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ + /*for (i=0; i<=cptcode; i++) {*/ + for (j=modmincovj; j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */ + printf("Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); + fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]); + if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */ + if( j != -1){ + ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th + covariate for which somebody answered excluding + undefined. Usually 2: 0 and 1. */ + } + ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th + covariate for which somebody answered including + undefined. Usually 3: -1, 0 and 1. */ + } /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for + * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ + } /* Ndum[-1] number of undefined modalities */ + + /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ + /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */ + /* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */ + /* modmincovj=3; modmaxcovj = 7; */ + /* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */ + /* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */ + /* defining two dummy variables: variables V1_1 and V1_2.*/ + /* nbcode[Tvar[j]][ij]=k; */ + /* nbcode[Tvar[j]][1]=0; */ + /* nbcode[Tvar[j]][2]=1; */ + /* 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*/ + 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*/ + cptcode = ij; /* New max modality for covar j */ + } /* end of loop on modality i=-1 to 1 or more */ + break; case 1: /* Testing on varying covariate, could be simple and * should look at waves or product of fixed * * varying. No time to test -1, assuming 0 and 1 only */ - ij=0; - for(i=0; i<=1;i++){ - nbcode[Tvar[k]][++ij]=i; - } - break; + ij=0; + for(i=0; i<=1;i++){ + nbcode[Tvar[k]][++ij]=i; + } + break; default: - break; + break; } /* end switch */ } /* end dummy test */ @@ -4733,22 +4773,22 @@ void concatwav(int wav[], int **dh, int ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */ Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/ Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */ - TmodelInvind[k]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ + TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ if(Fixed[k]!=0) anyvaryingduminmodel=1; - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ - /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ - /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ - /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ + /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ + /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ + /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ } } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ /* ij--; */ /* cptcoveff=ij; /\*Number of total covariates*\/ */ *cptcov=ij; /*Number of total real effective covariates: effective - * because they can be excluded from the model and real - * if in the model but excluded because missing values, but how to get k from ij?*/ + * because they can be excluded from the model and real + * if in the model but excluded because missing values, but how to get k from ij?*/ for(j=ij+1; j<= cptcovt; j++){ Tvaraff[j]=0; Tmodelind[j]=0; @@ -7876,21 +7916,21 @@ int decodemodel( char model[], int lasto if ((strpt=strstr(model,"age*age")) !=0){ printf(" strpt=%s, model=%s\n",strpt, model); if(strpt != model){ - printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ corresponding column of parameters.\n",model); - fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ corresponding column of parameters.\n",model); fflush(ficlog); - return 1; + return 1; } nagesqr=1; if (strstr(model,"+age*age") !=0) - substrchaine(modelsav, model, "+age*age"); + substrchaine(modelsav, model, "+age*age"); else if (strstr(model,"age*age+") !=0) - substrchaine(modelsav, model, "age*age+"); + substrchaine(modelsav, model, "age*age+"); else - substrchaine(modelsav, model, "age*age"); + substrchaine(modelsav, model, "age*age"); }else nagesqr=0; if (strlen(modelsav) >1){ @@ -8061,139 +8101,185 @@ Fixed[k] 0=fixed (product or simple), 1 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, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ - if (Tvar[k] <=ncovcol && (Typevar[k]==0 || Typevar[k]==2)){ /* Simple or product fixed dummy covariatee */ + if (Tvar[k] <=ncovcol && (Typevar[k]==0 || Typevar[k]==2)){ /* Simple or product fixed dummy (<=ncovcol) covariates */ Fixed[k]= 0; Dummy[k]= 0; ncoveff++; + modell[k].maintype= FTYPE; TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ - TvarFDind[ncoveff]=Tvar[k]; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ - }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ + TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */ Fixed[k]= 0; Dummy[k]= 1; nqfveff++; - TvarFQ[nqfveff]=Tvar[k]; /* TvarFQ[1]=V2 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ + modell[k].maintype= FTYPE; + modell[k].subtype= FQ; + TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){ Fixed[k]= 1; Dummy[k]= 0; ntveff++; /* Only simple time varying dummy variable */ - TvarVD[ntvveff]=Tvar[k]; /* TvarVD[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ - TvarVDind[ntveff++]=k; /* TvarVDind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ + modell[k].maintype= VTYPE; + modell[k].subtype= VD; + TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ + TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv); printf("Quasi TmodelInvind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv); - }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){ - Fixed[k]= 1; - Dummy[k]= 1; - TmodelInvQind[++nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */ - /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */ - printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv); + }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){ /* Only simple time varying quantitative variable V5*/ + Fixed[k]= 1; + Dummy[k]= 1; + nqtveff++; + modell[k].maintype= VTYPE; + modell[k].subtype= VQ; + TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ + TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ + TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */ + /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */ + printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv); printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv); }else if (Typevar[k] == 1) { /* product with age */ - if (Tvar[k] <=ncovcol ){ /* Simple or product fixed dummy covariatee */ - Fixed[k]= 2; - Dummy[k]= 2; - /* ncoveff++; */ + if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */ + Fixed[k]= 2; + Dummy[k]= 2; + modell[k].maintype= ATYPE; + modell[k].subtype= APFD; + /* ncoveff++; */ }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/ - Fixed[k]= 2; - Dummy[k]= 3; - /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */ + Fixed[k]= 2; + Dummy[k]= 3; + modell[k].maintype= ATYPE; + modell[k].subtype= APFQ; /* Product age * fixed quantitative */ + /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */ }else if( Tvar[k] <=ncovcol+nqv+ntv ){ - Fixed[k]= 3; - Dummy[k]= 2; - /* ntveff++; /\* Only simple time varying dummy variable *\/ */ + Fixed[k]= 3; + Dummy[k]= 2; + modell[k].maintype= ATYPE; + modell[k].subtype= APVD; /* Product age * varying dummy */ + /* ntveff++; /\* Only simple time varying dummy variable *\/ */ }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 3; - Dummy[k]= 3; - /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */ + Fixed[k]= 3; + Dummy[k]= 3; + modell[k].maintype= ATYPE; + modell[k].subtype= APVQ; /* Product age * varying quantitative */ + /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */ } }else if (Typevar[k] == 2) { /* product without age */ k1=Tposprod[k]; if(Tvard[k1][1] <=ncovcol){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 1; - Dummy[k]= 0; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[k]= 0; /* or 2 ?*/ - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 0; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 0; /* or 2 ?*/ + Dummy[k]= 1; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */ + } }else if(Tvard[k1][1] <=ncovcol+nqv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 0; /* or 2 ?*/ - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[k]= 0; /* or 2 ?*/ - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 0; /* or 2 ?*/ + Dummy[k]= 1; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */ + } }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 1; - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[k]= 1; - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 0; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 0; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */ + } }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[k]= 1; - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[k]= 1; - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[k]= 1; - Dummy[k]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[k]= 1; - Dummy[k]= 1; - } + if(Tvard[k1][2] <=ncovcol){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */ + }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ + Fixed[k]= 1; + Dummy[k]= 1; + modell[k].maintype= VTYPE; + modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */ + } }else{ - printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); - fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); + printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); + fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); } /* end k1 */ }else{ printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); } printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]); + printf(" modell[%d].maintype=%d, modell[%d].subtype=%d\n",k,modell[k].maintype,k,modell[k].subtype); fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]); } /* Searching for doublons in the model */ for(k1=1; k1<= cptcovt;k1++){ for(k2=1; k2