--- imach/src/imach.c 2016/07/12 18:42:34 1.226 +++ imach/src/imach.c 2016/08/23 07:40:50 1.233 @@ -1,6 +1,24 @@ -/* $Id: imach.c,v 1.226 2016/07/12 18:42:34 brouard Exp $ +/* $Id: imach.c,v 1.233 2016/08/23 07:40:50 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.233 2016/08/23 07:40:50 brouard + Summary: not working + + Revision 1.232 2016/08/22 14:20:21 brouard + Summary: not working + + 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 + + Revision 1.229 2016/07/23 09:45:53 brouard + Summary: Completing for func too + + Revision 1.228 2016/07/22 17:45:30 brouard + Summary: Fixing some arrays, still debugging + Revision 1.226 2016/07/12 18:42:34 brouard Summary: temp @@ -639,24 +657,27 @@ Short summary of the programme: - This program computes Healthy Life Expectancies from - cross-longitudinal data. Cross-longitudinal data consist in: -1- a - first survey ("cross") where individuals from different ages are - interviewed on their health status or degree of disability (in the - case of a health survey which is our main interest) -2- at least a - second wave of interviews ("longitudinal") which measure each change - (if any) in individual health status. Health expectancies are - computed from the time spent in each health state according to a - model. More health states you consider, more time is necessary to reach the - Maximum Likelihood of the parameters involved in the model. The - simplest model is the multinomial logistic model where pij is the - probability to be observed in state j at the second wave - conditional to be observed in state i at the first wave. Therefore - the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where - 'age' is age and 'sex' is a covariate. If you want to have a more - complex model than "constant and age", you should modify the program - where the markup *Covariates have to be included here again* invites - you to do it. More covariates you add, slower the + This program computes Healthy Life Expectancies or State-specific + (if states aren't health statuses) Expectancies from + cross-longitudinal data. Cross-longitudinal data consist in: + + -1- a first survey ("cross") where individuals from different ages + are interviewed on their health status or degree of disability (in + the case of a health survey which is our main interest) + + -2- at least a second wave of interviews ("longitudinal") which + measure each change (if any) in individual health status. Health + expectancies are computed from the time spent in each health state + according to a model. More health states you consider, more time is + necessary to reach the Maximum Likelihood of the parameters involved + in the model. The simplest model is the multinomial logistic model + where pij is the probability to be observed in state j at the second + wave conditional to be observed in state i at the first + wave. Therefore the model is: log(pij/pii)= aij + bij*age+ cij*sex + + etc , where 'age' is age and 'sex' is a covariate. If you want to + have a more complex model than "constant and age", you should modify + the program where the markup *Covariates have to be included here + again* invites you to do it. More covariates you add, slower the convergence. The advantage of this computer programme, compared to a simple @@ -678,21 +699,35 @@ of the life expectancies. It also computes the period (stable) prevalence. Back prevalence and projections: - - back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj) - Computes the back prevalence limit for any combination of covariate values k - at any age between ageminpar and agemaxpar and returns it in **bprlim. In the loops, - - **bprevalim(**bprlim, ***mobaverage, nlstate, *p, age, **oldm, **savm, **dnewm, **doldm, **dsavm, ftolpl, ncvyearp, k); - - hBijx Back Probability to be in state i at age x-h being in j at x + + - back_prevalence_limit(double *p, double **bprlim, double ageminpar, + double agemaxpar, double ftolpl, int *ncvyearp, double + dateprev1,double dateprev2, int firstpass, int lastpass, int + mobilavproj) + + Computes the back prevalence limit for any combination of + covariate values k at any age between ageminpar and agemaxpar and + returns it in **bprlim. In the loops, + + - **bprevalim(**bprlim, ***mobaverage, nlstate, *p, age, **oldm, + **savm, **dnewm, **doldm, **dsavm, ftolpl, ncvyearp, k); + + - hBijx Back Probability to be in state i at age x-h being in j at x Computes for any combination of covariates k and any age between bage and fage p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); oldm=oldms;savm=savms; - - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); + + - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' months (i.e. until age (in years) age+nhstepm*hstepm*stepm/12) by multiplying - nhstepm*hstepm matrices. Returns p3mat[i][j][h] after calling - p3mat[i][j][h]=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ - 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); + nhstepm*hstepm matrices. + + Returns p3mat[i][j][h] after calling + p3mat[i][j][h]=matprod2(newm, + bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, + dsavm,ij),\ 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, + oldm); Important routines @@ -864,12 +899,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.226 2016/07/12 18:42:34 brouard Exp $ */ +/* $Id: imach.c,v 1.233 2016/08/23 07:40:50 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.226 $ $Date: 2016/07/12 18:42:34 $"; +char fullversion[]="$Revision: 1.233 $ $Date: 2016/08/23 07:40:50 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -882,7 +917,11 @@ int cptcovsnq=0; /**< cptcovsnq number o int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */ int cptcovprodnoage=0; /**< Number of covariate products without age */ int cptcoveff=0; /* Total number of covariates to vary for printing results */ -int ncoveff=0; /* Total number of effective covariates in the model */ +int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */ +int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */ +int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */ + +int ncoveff=0; /* Total number of effective fixed dummy covariates in the model */ int nqfveff=0; /**< nqfveff Number of Quantitative Fixed Variables Effective */ int ntveff=0; /**< ntveff number of effective time varying variables */ int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */ @@ -907,6 +946,8 @@ int **dh; /* dh[mi][i] is number of step int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between * wave mi and wave mi+1 is not an exact multiple of stepm. */ int countcallfunc=0; /* Count the number of calls to func */ +int selected(int kvar); /* Is covariate kvar selected for printing results */ + double jmean=1; /* Mean space between 2 waves */ double **matprod2(); /* test */ double **oldm, **newm, **savm; /* Working pointers to matrices */ @@ -1036,15 +1077,70 @@ double ***cotvar; /* Time varying covari double ***cotqvar; /* Time varying quantitative covariate itqv */ double idx; int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ +int *TvarF; /**< TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +int *TvarFind; /**< TvarFind[1]=6, TvarFind[2]=7, Tvarind[3]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +int *TvarV; /**< TvarV[1]=Tvar[1]=5, TvarV[2]=Tvar[2]=4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +int *TvarVind; /**< TvarVind[1]=1, TvarVind[2]=2 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +int *TvarA; /**< TvarA[1]=Tvar[5]=5, TvarA[2]=Tvar[8]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +int *TvarAind; /**< TvarindA[1]=5, TvarAind[2]=8 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +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 */ -int *Fixed; /** Fixed[Tvar[k]] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */ -int *Dummy; /** Dummy[Tvar[k]] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */ +int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */ +int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */ int *Tage; +int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */ +int *Tmodelind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ +int *TmodelInvind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ +int *TmodelInvQind; /** Tmodelqind[1]=1 for V5(quantitative varying) position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ int *Ndum; /** Freq of modality (tricode */ /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ -int **Tvard, *Tprod, cptcovprod, *Tvaraff, *invalidvarcomb; +int **Tvard; +int *Tprod;/**< Gives the k position of the k1 product */ +int *Tposprod; /**< Gives the k1 product from the k position */ +/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3 + if V2+V1+V1*V4+age*V3+V3*V2 TProd[k1=2]=5 (V3*V2) + Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2 +*/ +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 */ @@ -2754,19 +2850,19 @@ double ***hpxij(double ***po, int nhstep agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ cov[2]=agexact; if(nagesqr==1) - cov[3]= agexact*agexact; + cov[3]= agexact*agexact; for (k=1; k<=cptcovn;k++) - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; - /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ - /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; - /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ + /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; - /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ - - + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ + + /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ /* right multiplication of oldm by the current matrix */ @@ -2917,10 +3013,9 @@ double func( double *x) int ioffset=0; double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; double **out; - double sw; /* Sum of weights */ double lli; /* Individual log likelihood */ int s1, s2; - int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quatitative time varying covariate */ + int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quantitative time varying covariate */ double bbh, survp; long ipmx; double agexact; @@ -2945,14 +3040,10 @@ double func( double *x) to be observed in j being in i according to the model. */ ioffset=2+nagesqr+cptcovage; - /* for (k=1; k<=cptcovn;k++){ /\* Simple and product covariates without age* products *\/ */ - for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */ - cov[++ioffset]=covar[Tvar[k]][i]; - } - for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitatives and Fixed covariates */ - cov[++ioffset]=coqvar[iqv][i]; - } - + /* Fixed */ + for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */ + 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)*/ + } /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] has been calculated etc */ @@ -2966,92 +3057,85 @@ 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]][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]][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 */ @@ -3246,15 +3330,16 @@ double func( double *x) /*************** log-likelihood *************/ double funcone( double *x) { - /* Same as likeli but slower because of a lot of printf and if */ + /* Same as func but slower because of a lot of printf and if */ int i, ii, j, k, mi, d, kk; - int ioffset=0; + int ioffset=0; double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; double **out; double lli; /* Individual log likelihood */ double llt; int s1, s2; - int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate */ + int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quantitative time varying covariate */ + double bbh, survp; double agexact; double agebegin, ageend; @@ -3270,48 +3355,72 @@ double funcone( double *x) ioffset=0; for (i=1,ipmx=0, sw=0.; i<=imx; i++){ ioffset=2+nagesqr+cptcovage; + /* Fixed */ /* 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 covariates without age* products */ - cov[++ioffset]=covar[Tvar[k]][i]; - } - for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitative fixed covariates */ - cov[++ioffset]=coqvar[Tvar[iqv]][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 */ + 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]; */ +/* cov[2+6]=covar[2][i]; V2 */ +/* cov[TvarFind[2]]=covar[Tvar[TvarFind[2]]][i]; */ +/* cov[2+7]=covar[Tvar[7]][i]; */ +/* cov[2+7]=covar[7][i]; V7=V1*V2 */ +/* cov[TvarFind[3]]=covar[Tvar[TvarFind[3]]][i]; */ +/* cov[2+9]=covar[Tvar[9]][i]; */ +/* cov[2+9]=covar[1][i]; V1 */ } + /* 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 */ - cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i]; - } - for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */ - cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i]; - } + /* Wave varying (but not age varying) */ + for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ + cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][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]; */ + /* } */ for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ for(d=0; d nlstate && (mle <5) ){ /* Jackson */ - lli=log(out[s1][s2] - savm[s1][s2]); + lli=log(out[s1][s2] - savm[s1][s2]); } else if ( s2==-1 ) { /* alive */ - for (j=1,survp=0. ; j<=nlstate; j++) - survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; - lli= log(survp); + for (j=1,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); }else if (mle==1){ - lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ } else if(mle==2){ - lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ + lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ } else if(mle==3){ /* exponential inter-extrapolation */ - lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ + lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ } else if (mle==4){ /* mle=4 no inter-extrapolation */ - lli=log(out[s1][s2]); /* Original formula */ + lli=log(out[s1][s2]); /* Original formula */ } else{ /* mle=0 back to 1 */ - lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ - /*lli=log(out[s1][s2]); */ /* Original formula */ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + /*lli=log(out[s1][s2]); */ /* Original formula */ } /* End of if */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ if(globpr){ - fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ + fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ - num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, - 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); - for(k=1,llt=0.,l=0.; k<=nlstate; k++){ - llt +=ll[k]*gipmx/gsw; - fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); - } - fprintf(ficresilk," %10.6f\n", -llt); + num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, + 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); + for(k=1,llt=0.,l=0.; k<=nlstate; k++){ + llt +=ll[k]*gipmx/gsw; + fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); + } + fprintf(ficresilk," %10.6f\n", -llt); } - } /* end of wave */ - } /* end of individual */ - for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; - /* printf("l1=%f l2=%f ",ll[1],ll[2]); */ - l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ - if(globpr==0){ /* First time we count the contributions and weights */ - gipmx=ipmx; - gsw=sw; - } - return -l; + } /* end of wave */ +} /* end of individual */ +for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; +/* printf("l1=%f l2=%f ",ll[1],ll[2]); */ +l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ +if(globpr==0){ /* First time we count the contributions and weights */ + gipmx=ipmx; + gsw=sw; +} +return -l; } @@ -3867,7 +3976,7 @@ void freqsummary(char fileres[], int ia int firstpass, int lastpass, int stepm, int weightopt, char model[]) { /* Some frequencies */ - int i, m, jk, j1, bool, z1,j; + int i, m, jk, j1, bool, z1,j, k, iv; int iind=0, iage=0; int mi; /* Effective wave */ int first; @@ -3928,7 +4037,8 @@ Title=%s
Datafile=%s Firstpass=%d La freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); j1=0; - j=ncoveff; + /* j=ncoveff; /\* Only fixed dummy covariates *\/ */ + j=cptcoveff; /* Only dummy covariates of the model */ if (cptcovn<1) {j=1;ncodemax[1]=1;} first=1; @@ -3940,108 +4050,125 @@ Title=%s
Datafile=%s Firstpass=%d La Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff */ - for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination excluding varying and quantitatives */ + for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */ posproptt=0.; /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); 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; } - for (z1=1; z1<= nqfveff; z1++) { - meanq[z1]+=0.; - for(m=1;m<=lastpass;m++){ - meanqt[m][z1]=0.; - } - } - + /* for (z1=1; z1<= nqfveff; z1++) { */ + /* meanq[z1]+=0.; */ + /* for(m=1;m<=lastpass;m++){ */ + /* meanqt[m][z1]=0.; */ + /* } */ + /* } */ + dateintsum=0; k2cpt=0; - /* For that comination of covariate j1, we count and print the frequencies */ + /* For 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 (nqfveff >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]; - } - for (z1=1; z1<=ncoveff; z1++) { - /* if(Tvaraff[z1] ==-20){ */ - /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ - /* }else if(Tvaraff[z1] ==-10){ */ - /* /\* sumnew+=coqvar[z1][iind]; *\/ */ - /* }else */ - if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ - /* Tests if this individual i responded to j1 (V4=1 V3=0) */ - bool=0; - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", - bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), - 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*/ - } - } /* end z1 */ - } /* cptcovn > 0 */ - - if (bool==1){ /* We selected an individual iin satisfying combination j1 */ - /* 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 m */ + 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++) { + /* 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 */ } /* 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) { + /* if (ncoveff>0) { */ + if (cptcoveff>0) { fprintf(ficresp, "\n#********** Variable "); fprintf(ficresphtm, "\n

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

********** Variable "); - for (z1=1; z1<=ncoveff; z1++){ - fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + 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, "**********\n#"); fprintf(ficresphtm, "**********

\n"); fprintf(ficresphtmfr, "**********\n"); fprintf(ficlog, "\n#********** Variable "); - for (z1=1; z1<=ncoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficlog, "**********\n"); } fprintf(ficresphtm,""); @@ -4051,126 +4178,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"); @@ -4188,7 +4315,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); @@ -4203,97 +4330,107 @@ Title=%s
Datafile=%s Firstpass=%d La } /************ Prevalence ********************/ - void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) - { - /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people - in each health status at the date of interview (if between dateprev1 and dateprev2). - We still use firstpass and lastpass as another selection. - */ +void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) +{ + /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people + in each health status at the date of interview (if between dateprev1 and dateprev2). + We still use firstpass and lastpass as another selection. + */ - int i, m, jk, j1, bool, z1,j; - int mi; /* Effective wave */ - int iage; - double agebegin, ageend; - - double **prop; - double posprop; - double y2; /* in fractional years */ - int iagemin, iagemax; - int first; /** to stop verbosity which is redirected to log file */ - - iagemin= (int) agemin; - iagemax= (int) agemax; - /*pp=vector(1,nlstate);*/ - prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); - /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ - j1=0; + int i, m, jk, j1, bool, z1,j, iv; + int mi; /* Effective wave */ + int iage; + double agebegin, ageend; + + double **prop; + double posprop; + double y2; /* in fractional years */ + int iagemin, iagemax; + int first; /** to stop verbosity which is redirected to log file */ + + iagemin= (int) agemin; + iagemax= (int) agemax; + /*pp=vector(1,nlstate);*/ + prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); + /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ + j1=0; - /*j=cptcoveff;*/ - if (cptcovn<1) {j=1;ncodemax[1]=1;} + /*j=cptcoveff;*/ + if (cptcovn<1) {j=1;ncodemax[1]=1;} - first=1; - for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ - for (i=1; i<=nlstate; i++) - for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) - prop[i][iage]=0.0; - - for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; - if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/ - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) - bool=0; - } - if (bool==1) { /* For this combination of covariates values, this individual fits */ - /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ - for(mi=1; mi=firstpass && m <=lastpass){ - y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ - if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ - if(agev[m][i]==0) agev[m][i]=iagemax+1; - if(agev[m][i]==1) agev[m][i]=iagemax+2; - if((int)agev[m][i] iagemax+3+AGEMARGE){ - printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); - exit(1); - } - if (s[m][i]>0 && s[m][i]<=nlstate) { - /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ - prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ - prop[s[m][i]][iagemax+3] += weight[i]; - } /* end valid statuses */ - } /* end selection of dates */ - } /* end selection of waves */ - } /* end effective waves */ - } /* end bool */ - } - for(i=iagemin; i <= iagemax+3; i++){ - for(jk=1,posprop=0; jk <=nlstate ; jk++) { - posprop += prop[jk][i]; - } - - for(jk=1; jk <=nlstate ; jk++){ - if( i <= iagemax){ - if(posprop>=1.e-5){ - probs[i][jk][j1]= prop[jk][i]/posprop; - } else{ - if(first==1){ - first=0; - printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); - } - } - } - }/* end jk */ - }/* end i */ + first=1; + 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+3+AGEMARGE; iage++) + prop[i][iage]=0.0; + printf("Prevalence combination of varying and fixed dummies %d\n",j1); + /* fprintf(ficlog," V%d=%d ",Tvaraff[j1],nbcode[Tvaraff[j1]][codtabm(k,j1)]); */ + fprintf(ficlog,"Prevalence combination of varying and fixed dummies %d\n",j1); + + for (i=1; i<=imx; i++) { /* Each individual */ + bool=1; + /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ + for(mi=1; mi=firstpass && m <=lastpass){ + y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ + if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ + if(agev[m][i]==0) agev[m][i]=iagemax+1; + if(agev[m][i]==1) agev[m][i]=iagemax+2; + if((int)agev[m][i] iagemax+3+AGEMARGE){ + printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); + exit(1); + } + if (s[m][i]>0 && s[m][i]<=nlstate) { + /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ + prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ + prop[s[m][i]][iagemax+3] += weight[i]; + } /* end valid statuses */ + } /* end selection of dates */ + } /* end selection of waves */ + } /* end bool */ + } /* end wave */ + } /* end individual */ + for(i=iagemin; i <= iagemax+3; i++){ + for(jk=1,posprop=0; jk <=nlstate ; jk++) { + posprop += prop[jk][i]; + } + + for(jk=1; jk <=nlstate ; jk++){ + if( i <= iagemax){ + if(posprop>=1.e-5){ + probs[i][jk][j1]= prop[jk][i]/posprop; + } else{ + if(first==1){ + first=0; + printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,j1,probs[i][jk][j1]); + } + } + } + }/* end jk */ + }/* end i */ /*} *//* end i1 */ - } /* end j1 */ + } /* end j1 */ - /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ - /*free_vector(pp,1,nlstate);*/ - free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE); - } /* End of prevalence */ + /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ + /*free_vector(pp,1,nlstate);*/ + free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE); +} /* End of prevalence */ /************* Waves Concatenation ***************/ @@ -4304,7 +4441,7 @@ void concatwav(int wav[], int **dh, int mw[mi][i] is the mi (mi=1 to wav[i]) effective wave of individual i dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] and mw[mi+1][i]. dh depends on stepm. - */ + */ int i=0, mi=0, m=0, mli=0; /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; @@ -4323,85 +4460,85 @@ void concatwav(int wav[], int **dh, int /* Treating live states */ for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ mi=0; /* First valid wave */ - mli=0; /* Last valid wave */ + mli=0; /* Last valid wave */ m=firstpass; while(s[m][i] <= nlstate){ /* a live 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; - mli=m; + if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */ + mli=m-1;/* mw[++mi][i]=m-1; */ + }else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ + mw[++mi][i]=m; + mli=m; } /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */ if(m < lastpass){ /* m < lastpass, standard case */ - m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */ + 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; + break; #else - if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ - 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 pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); - firsthree=1; - } - fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); - mw[++mi][i]=m; - mli=m; - } - if(s[m][i]==-2){ /* Vital status is really unknown */ - nbwarn++; - if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ - printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); - fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); - } - break; - } - break; + if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ + if(firsthree == 0){ + printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + firsthree=1; + } + fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + mw[++mi][i]=m; + mli=m; + } + if(s[m][i]==-2){ /* Vital status is really unknown */ + nbwarn++; + if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ + printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); + fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); + } + break; + } + break; #endif - }/* End m >= lastpass */ + }/* End m >= lastpass */ }/* end while */ - /* mi is the last effective wave, m is lastpass, mw[j][i] gives the # of j-th effective wave for individual i */ + /* mi is the last effective wave, m is lastpass, mw[j][i] gives the # of j-th effective wave for individual i */ /* After last pass */ /* Treating death states */ if (s[m][i] > nlstate){ /* In a death state */ - /* if( mint[m][i]==mdc[m][i] && anint[m][i]==andc[m][i]){ /\* same date of death and date of interview *\/ */ - /* } */ + /* if( mint[m][i]==mdc[m][i] && anint[m][i]==andc[m][i]){ /\* same date of death and date of interview *\/ */ + /* } */ mi++; /* Death is another wave */ /* if(mi==0) never been interviewed correctly before death */ - /* Only death is a correct wave */ + /* Only death is a correct wave */ mw[mi][i]=m; } #ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE - else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */ + else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */ /* m++; */ /* mi++; */ /* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */ /* mw[mi][i]=m; */ 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 */ - nbwarn++; - if(firstfiv==0){ - printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); - firstfiv=1; - }else{ - fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); - } - }else{ /* Death occured afer last wave potential bias */ - nberr++; - if(firstwo==0){ - printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); - firstwo=1; - } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); - } + if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */ + nbwarn++; + if(firstfiv==0){ + printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + firstfiv=1; + }else{ + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + } + }else{ /* Death occured afer last wave potential bias */ + nberr++; + if(firstwo==0){ + printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + firstwo=1; + } + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + } }else{ /* end date of interview is known */ - /* death is known but not confirmed by death status at any wave */ - if(firstfour==0){ - printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); - firstfour=1; - } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + /* death is known but not confirmed by death status at any wave */ + if(firstfour==0){ + printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + firstfour=1; + } + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); } } /* end if date of death is known */ #endif @@ -4410,11 +4547,11 @@ void concatwav(int wav[], int **dh, int if(mi==0){ nbwarn++; if(first==0){ - printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); - first=1; + printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); + first=1; } if(first==1){ - fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); + fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); } } /* end mi==0 */ } /* End individuals */ @@ -4424,92 +4561,92 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=imx; i++){ for(mi=1; mi nlstate) { /* A death */ - if (agedc[i] < 2*AGESUP) { - j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); - if(j==0) j=1; /* Survives at least one month after exam */ - else if(j<0){ - nberr++; - printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - j=1; /* Temporary Dangerous patch */ - printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); - fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); - } - k=k+1; - if (j >= jmax){ - jmax=j; - ijmax=i; - } - if (j <= jmin){ - jmin=j; - ijmin=i; - } - sum=sum+j; - /*if (j<0) printf("j=%d num=%d \n",j,i);*/ - /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ - } - } - else{ - j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); + if (s[mw[mi+1][i]][i] > nlstate) { /* A death */ + if (agedc[i] < 2*AGESUP) { + j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); + if(j==0) j=1; /* Survives at least one month after exam */ + else if(j<0){ + nberr++; + printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + j=1; /* Temporary Dangerous patch */ + printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); + fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); + } + k=k+1; + if (j >= jmax){ + jmax=j; + ijmax=i; + } + if (j <= jmin){ + jmin=j; + ijmin=i; + } + sum=sum+j; + /*if (j<0) printf("j=%d num=%d \n",j,i);*/ + /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ + } + } + else{ + j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); /* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ - k=k+1; - if (j >= jmax) { - jmax=j; - ijmax=i; - } - else if (j <= jmin){ - jmin=j; - ijmin=i; - } - /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ - /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ - if(j<0){ - nberr++; - printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - } - sum=sum+j; - } - jk= j/stepm; - jl= j -jk*stepm; - ju= j -(jk+1)*stepm; - if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */ - if(jl==0){ - dh[mi][i]=jk; - bh[mi][i]=0; - }else{ /* We want a negative bias in order to only have interpolation ie - * to avoid the price of an extra matrix product in likelihood */ - dh[mi][i]=jk+1; - bh[mi][i]=ju; - } - }else{ - if(jl <= -ju){ - dh[mi][i]=jk; - bh[mi][i]=jl; /* bias is positive if real duration - * is higher than the multiple of stepm and negative otherwise. - */ - } - else{ - dh[mi][i]=jk+1; - bh[mi][i]=ju; - } - if(dh[mi][i]==0){ - dh[mi][i]=1; /* At least one step */ - bh[mi][i]=ju; /* At least one step */ - /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/ - } - } /* end if mle */ + k=k+1; + if (j >= jmax) { + jmax=j; + ijmax=i; + } + else if (j <= jmin){ + jmin=j; + ijmin=i; + } + /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ + /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ + if(j<0){ + nberr++; + printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + } + sum=sum+j; + } + jk= j/stepm; + jl= j -jk*stepm; + ju= j -(jk+1)*stepm; + if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */ + if(jl==0){ + dh[mi][i]=jk; + bh[mi][i]=0; + }else{ /* We want a negative bias in order to only have interpolation ie + * to avoid the price of an extra matrix product in likelihood */ + dh[mi][i]=jk+1; + bh[mi][i]=ju; + } + }else{ + if(jl <= -ju){ + dh[mi][i]=jk; + bh[mi][i]=jl; /* bias is positive if real duration + * is higher than the multiple of stepm and negative otherwise. + */ + } + else{ + dh[mi][i]=jk+1; + bh[mi][i]=ju; + } + if(dh[mi][i]==0){ + dh[mi][i]=1; /* At least one step */ + bh[mi][i]=ju; /* At least one step */ + /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/ + } + } /* end if mle */ } } /* end wave */ } jmean=sum/k; printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean); fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %d) Max=%d (%d) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean); - } +} /*********** Tricode ****************************/ void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum) @@ -4534,82 +4671,90 @@ void concatwav(int wav[], int **dh, int /* Loop on covariates without age and products and no quantitative variable */ /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */ - for (j=1; j<=cptcovsnq; j++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ - for (k=-1; k < maxncov; k++) Ndum[k]=0; - for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the - modality of this covariate Vj*/ - switch(Typevar[j]) { - case 1: /* A real fixed dummy covariate */ - ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i - * If product of Vn*Vm, still boolean *: - * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables - * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */ - /* Finds for covariate j, n=Tvar[j] of Vn . ij is the - modality of the nth covariate of individual i. */ - if (ij > modmaxcovj) - modmaxcovj=ij; - else if (ij < modmincovj) - modmincovj=ij; - if ((ij < -1) && (ij > NCOVMAX)){ - printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); - exit(1); - }else - Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/ - /* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */ - /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ - /* getting the maximum value of the modality of the covariate - (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and - female ies 1, then modmaxcovj=1.*/ - break; - case 2: - break; - - } - } /* end for loop on individuals i */ - printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); - fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); - cptcode=modmaxcovj; - /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ - /*for (i=0; i<=cptcode; i++) {*/ - for (k=modmincovj; k<=modmaxcovj; k++) { /* k=-1 ? 0 and 1*//* For each value k of the modality of model-cov j */ - printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); - fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); - if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */ - if( k != -1){ - ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th - covariate for which somebody answered excluding - undefined. Usually 2: 0 and 1. */ - } - ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th - covariate for which somebody answered including - undefined. Usually 3: -1, 0 and 1. */ - } - /* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for - * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ - } /* 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[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ - cptcode = ij; /* New max modality for covar j */ - } /* end of loop on modality i=-1 to 1 or more */ + for (k=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 */ + 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; + 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; + default: + 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 *\/ */ @@ -4627,32 +4772,50 @@ void concatwav(int wav[], int **dh, int } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; - + /* Look at fixed dummy (single or product) covariates to check empty modalities */ for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ - ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ - Ndum[ij]++; /* Might be supersed V1 + V1*age */ + ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */ + Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */ + /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, {2, 1, 1, 1, 2, 1, 1, 0, 0} */ } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */ ij=0; - for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ + /* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */ + for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ - if((Ndum[i]!=0) && (i<=ncovcol)){ + /* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */ + if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */ + /* If product not in single variable we don't print results */ /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ - Tvaraff[++ij]=i; /*For printing (unclear) */ - }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ - Tvaraff[++ij]=-10; /* Dont'n know how to treat quantitative variables yet */ - }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ - Tvaraff[++ij]=i; /*For printing (unclear) */ - }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ - Tvaraff[++ij]=-20; /* Dont'n know how to treat quantitative variables yet */ - } + ++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[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 *\/ */ + } } /* 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*/ + * 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; + } + for(j=ntveff+1; j<= cptcovt; j++){ + TmodelInvind[j]=0; + } + /* To be sorted */ + ; } @@ -5802,6 +5965,7 @@ void printinghtml(char fileresu[], char fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); } + /* if(nqfveff+nqtveff 0) */ /* Test to be done */ fprintf(fichtm," ************\n
"); if(invalidvarcomb[k1]){ fprintf(fichtm,"\n

Combination (%d) ignored because no cases

\n",k1); @@ -5997,7 +6161,7 @@ void printinggnuplot(char fileresu[], ch strcpy(optfileres,"vpl"); /* 1eme*/ for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ - for (k1=1; k1<= m ; k1 ++) { /* For each valid combination of covariate */ + for (k1=1; k1<= m && selected(k1) ; k1 ++) { /* For each valid combination of covariate */ /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ @@ -6056,8 +6220,8 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ if(k==cptcoveff){ - fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' with line ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ - 6+(cpt-1), cpt ); + fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ + 4+(cpt-1), cpt ); /* 4 or 6 ?*/ }else{ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]); kl++; @@ -6210,133 +6374,133 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ - lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ - /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ - /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ - /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ - vlv= nbcode[Tvaraff[k]][lv]; - fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[k]][lv]; + fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); } fprintf(ficgp,"\n#\n"); if(invalidvarcomb[k1]){ - fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); - continue; - } - + fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); + continue; + } + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ -set ter svg size 640, 480\n \ -unset log y\n \ +set ter svg size 640, 480\n \ +unset log y\n \ plot [%.f:%.f] ", ageminpar, agemaxpar); k=3; for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ - if(j==1) - fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); - else - fprintf(ficgp,", '' "); - l=(nlstate+ndeath)*(cpt-1) +j; - fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l); - /* for (i=2; i<= nlstate+ndeath ; i ++) */ - /* fprintf(ficgp,"+$%d",k+l+i-1); */ - fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j); + if(j==1) + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); + else + fprintf(ficgp,", '' "); + l=(nlstate+ndeath)*(cpt-1) +j; + fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l); + /* for (i=2; i<= nlstate+ndeath ; i ++) */ + /* fprintf(ficgp,"+$%d",k+l+i-1); */ + fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j); } /* nlstate */ fprintf(ficgp,", '' "); fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1); for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ - l=(nlstate+ndeath)*(cpt-1) +j; - if(j < nlstate) - fprintf(ficgp,"$%d +",k+l); - else - fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt); + l=(nlstate+ndeath)*(cpt-1) +j; + if(j < nlstate) + fprintf(ficgp,"$%d +",k+l); + else + fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt); } fprintf(ficgp,"\nset out\n"); } /* end cpt state*/ } /* end covariate */ - + /* 6eme */ /* CV preval stable (period) for each covariate */ for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ - + fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ - lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ - /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ - /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ - /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ - vlv= nbcode[Tvaraff[k]][lv]; - fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[k]][lv]; + fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); } fprintf(ficgp,"\n#\n"); if(invalidvarcomb[k1]){ - fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); - continue; + fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); + continue; } - + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ -set ter svg size 640, 480\n \ -unset log y\n \ +set ter svg size 640, 480\n \ +unset log y\n \ plot [%.f:%.f] ", ageminpar, agemaxpar); k=3; /* Offset */ for (i=1; i<= nlstate ; i ++){ - if(i==1) - fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); - else - fprintf(ficgp,", '' "); - l=(nlstate+ndeath)*(i-1)+1; - fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); - for (j=2; j<= nlstate ; j ++) - fprintf(ficgp,"+$%d",k+l+j-1); - fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt); + if(i==1) + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); + else + fprintf(ficgp,", '' "); + l=(nlstate+ndeath)*(i-1)+1; + fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); + for (j=2; j<= nlstate ; j ++) + fprintf(ficgp,"+$%d",k+l+j-1); + fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt); } /* nlstate */ fprintf(ficgp,"\nset out\n"); } /* end cpt state*/ } /* end covariate */ - - + + /* 7eme */ if(backcast == 1){ /* CV back preval stable (period) for each covariate */ for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ - fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); - for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ - lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ - /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ - /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); + for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ - vlv= nbcode[Tvaraff[k]][lv]; - fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); - } - fprintf(ficgp,"\n#\n"); - if(invalidvarcomb[k1]){ - fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); - continue; - } - - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ -set ter svg size 640, 480\n \ -unset log y\n \ + vlv= nbcode[Tvaraff[k]][lv]; + fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); + } + fprintf(ficgp,"\n#\n"); + if(invalidvarcomb[k1]){ + fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); + continue; + } + + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ +set ter svg size 640, 480\n \ +unset log y\n \ plot [%.f:%.f] ", ageminpar, agemaxpar); - k=3; /* Offset */ - for (i=1; i<= nlstate ; i ++){ - if(i==1) - fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_")); - else - fprintf(ficgp,", '' "); - /* l=(nlstate+ndeath)*(i-1)+1; */ - l=(nlstate+ndeath)*(cpt-1)+1; - /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */ - /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */ - fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */ - /* for (j=2; j<= nlstate ; j ++) */ - /* fprintf(ficgp,"+$%d",k+l+j-1); */ - /* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */ - fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt); - } /* nlstate */ - fprintf(ficgp,"\nset out\n"); + k=3; /* Offset */ + for (i=1; i<= nlstate ; i ++){ + if(i==1) + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_")); + else + fprintf(ficgp,", '' "); + /* l=(nlstate+ndeath)*(i-1)+1; */ + l=(nlstate+ndeath)*(cpt-1)+1; + /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */ + /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */ + fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */ + /* for (j=2; j<= nlstate ; j ++) */ + /* fprintf(ficgp,"+$%d",k+l+j-1); */ + /* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */ + fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt); + } /* nlstate */ + fprintf(ficgp,"\nset out\n"); } /* end cpt state*/ } /* end covariate */ } /* End if backcast */ @@ -6347,110 +6511,110 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ - fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); - for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ - lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ - /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ - /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ - /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ - vlv= nbcode[Tvaraff[k]][lv]; - fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); - } - fprintf(ficgp,"\n#\n"); - if(invalidvarcomb[k1]){ - fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); - continue; - } - - fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ -set ter svg size 640, 480\n \ -unset log y\n \ + fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); + for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[k]][lv]; + fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); + } + fprintf(ficgp,"\n#\n"); + if(invalidvarcomb[k1]){ + fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); + continue; + } + + fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ +set ter svg size 640, 480\n \ +unset log y\n \ plot [%.f:%.f] ", ageminpar, agemaxpar); - for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ - /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ - /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ - /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ - /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ - if(i==1){ - fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); - }else{ - fprintf(ficgp,",\\\n '' "); - } - if(cptcoveff ==0){ /* No covariate */ - ioffset=2; /* Age is in 2 */ - /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ - /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ - /*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ - /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ - fprintf(ficgp," u %d:(", ioffset); - if(i==nlstate+1) - fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ", \ - ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); - else - fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ - ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); - }else{ /* more than 2 covariates */ - if(cptcoveff ==1){ - ioffset=4; /* Age is in 4 */ - }else{ - ioffset=6; /* Age is in 6 */ - /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ - /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ - } - fprintf(ficgp," u %d:(",ioffset); - kl=0; - strcpy(gplotcondition,"("); - for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */ - lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */ - /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ - /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ - /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ - vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */ - kl++; - sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); - kl++; - if(k 1) - sprintf(gplotcondition+strlen(gplotcondition)," && "); - } - strcpy(gplotcondition+strlen(gplotcondition),")"); - /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ - /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ - /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ - /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ - if(i==nlstate+1){ - fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \ - ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); - }else{ - fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ - ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt ); - } - } /* end if covariate */ - } /* nlstate */ - fprintf(ficgp,"\nset out\n"); + for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ + /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ + /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ + if(i==1){ + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); + }else{ + fprintf(ficgp,",\\\n '' "); + } + if(cptcoveff ==0){ /* No covariate */ + ioffset=2; /* Age is in 2 */ + /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ + /*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ + fprintf(ficgp," u %d:(", ioffset); + if(i==nlstate+1) + fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ", \ + ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); + else + fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ + ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); + }else{ /* more than 2 covariates */ + if(cptcoveff ==1){ + ioffset=4; /* Age is in 4 */ + }else{ + ioffset=6; /* Age is in 6 */ + /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ + /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ + } + fprintf(ficgp," u %d:(",ioffset); + kl=0; + strcpy(gplotcondition,"("); + for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */ + lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */ + /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ + /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ + /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ + vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */ + kl++; + sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); + kl++; + if(k 1) + sprintf(gplotcondition+strlen(gplotcondition)," && "); + } + strcpy(gplotcondition+strlen(gplotcondition),")"); + /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ + /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ + /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ + /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ + if(i==nlstate+1){ + fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \ + ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); + }else{ + fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ + ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt ); + } + } /* end if covariate */ + } /* nlstate */ + fprintf(ficgp,"\nset out\n"); } /* end cpt state*/ } /* end covariate */ } /* End if prevfcast */ - - + + /* proba elementaires */ fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); for(i=1,jk=1; i <=nlstate; i++){ fprintf(ficgp,"# initial state %d\n",i); for(k=1; k <=(nlstate+ndeath); k++){ if (k != i) { - fprintf(ficgp,"# current state %d\n",k); - for(j=1; j <=ncovmodel; j++){ - fprintf(ficgp,"p%d=%f; ",jk,p[jk]); - jk++; - } - fprintf(ficgp,"\n"); + fprintf(ficgp,"# current state %d\n",k); + for(j=1; j <=ncovmodel; j++){ + fprintf(ficgp,"p%d=%f; ",jk,p[jk]); + jk++; + } + fprintf(ficgp,"\n"); } } } fprintf(ficgp,"##############\n#\n"); - + /*goto avoid;*/ fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n"); fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n"); @@ -6521,17 +6685,17 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) } } else - fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); /* Valgrind bug nbcode */ + fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); /* Valgrind bug nbcode */ } }else{ i=i-ncovmodel; if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */ fprintf(ficgp," (1."); } - + if(ng != 1){ fprintf(ficgp,")/(1"); - + for(k1=1; k1 <=nlstate; k1++){ if(nagesqr==0) fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); @@ -6784,66 +6948,67 @@ void prevforecast(char fileres[], double if(jprojmean==0) jprojmean=1; if(mprojmean==0) jprojmean=1; - i1=cptcoveff; + i1=pow(2,cptcoveff); if (cptcovn < 1){i1=1;} fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); fprintf(ficresf,"#****** Routine prevforecast **\n"); - + /* if (h==(int)(YEARM*yearp)){ */ - for(cptcov=1, k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ - k=k+1; - fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - } - fprintf(ficresf," yearproj age"); - for(j=1; j<=nlstate+ndeath;j++){ - for(i=1; i<=nlstate;i++) - fprintf(ficresf," p%d%d",i,j); - fprintf(ficresf," wp.%d",j); - } - for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { - fprintf(ficresf,"\n"); - fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); - for (agec=fage; agec>=(ageminpar-1); agec--){ - nhstepm=(int) rint((agelim-agec)*YEARM/stepm); - nhstepm = nhstepm/hstepm; - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - oldm=oldms;savm=savms; - hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k); - - for (h=0; h<=nhstepm; h++){ - if (h*hstepm/YEARM*stepm ==yearp) { - fprintf(ficresf,"\n"); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); - } - for(j=1; j<=nlstate+ndeath;j++) { - ppij=0.; - for(i=1; i<=nlstate;i++) { - if (mobilav==1) - ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; - else { - ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; - } - if (h*hstepm/YEARM*stepm== yearp) { - fprintf(ficresf," %.3f", p3mat[i][j][h]); - } - } /* end i */ - if (h*hstepm/YEARM*stepm==yearp) { - fprintf(ficresf," %.3f", ppij); - } - }/* end j */ - } /* end h */ - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - } /* end agec */ - } /* end yearp */ - } /* end cptcod */ - } /* end cptcov */ + for(k=1;k<=i1;k++){ + if(invalidvarcomb[k]){ + printf("\nCombination (%d) projection ignored because no cases \n",k); + continue; + } + fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + } + fprintf(ficresf," yearproj age"); + for(j=1; j<=nlstate+ndeath;j++){ + for(i=1; i<=nlstate;i++) + fprintf(ficresf," p%d%d",i,j); + fprintf(ficresf," wp.%d",j); + } + for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { + fprintf(ficresf,"\n"); + fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); + for (agec=fage; agec>=(ageminpar-1); agec--){ + nhstepm=(int) rint((agelim-agec)*YEARM/stepm); + nhstepm = nhstepm/hstepm; + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + oldm=oldms;savm=savms; + hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k); + + for (h=0; h<=nhstepm; h++){ + if (h*hstepm/YEARM*stepm ==yearp) { + fprintf(ficresf,"\n"); + for(j=1;j<=cptcoveff;j++) + fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); + } + for(j=1; j<=nlstate+ndeath;j++) { + ppij=0.; + for(i=1; i<=nlstate;i++) { + if (mobilav==1) + ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k]; + else { + ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; + } + if (h*hstepm/YEARM*stepm== yearp) { + fprintf(ficresf," %.3f", p3mat[i][j][h]); + } + } /* end i */ + if (h*hstepm/YEARM*stepm==yearp) { + fprintf(ficresf," %.3f", ppij); + } + }/* end j */ + } /* end h */ + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + } /* end agec */ + } /* end yearp */ + } /* end k */ fclose(ficresf); printf("End of Computing forecasting \n"); @@ -6982,165 +7147,165 @@ void prevforecast(char fileres[], double /* } */ /************** Forecasting *****not tested NB*************/ -void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ +/* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */ - int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; - int *popage; - double calagedatem, agelim, kk1, kk2; - double *popeffectif,*popcount; - double ***p3mat,***tabpop,***tabpopprev; - /* double ***mobaverage; */ - char filerespop[FILENAMELENGTH]; +/* int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; */ +/* int *popage; */ +/* double calagedatem, agelim, kk1, kk2; */ +/* double *popeffectif,*popcount; */ +/* double ***p3mat,***tabpop,***tabpopprev; */ +/* /\* double ***mobaverage; *\/ */ +/* char filerespop[FILENAMELENGTH]; */ - tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - agelim=AGESUP; - calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; +/* tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */ +/* tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */ +/* agelim=AGESUP; */ +/* calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; */ - prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); +/* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */ - strcpy(filerespop,"POP_"); - strcat(filerespop,fileresu); - if((ficrespop=fopen(filerespop,"w"))==NULL) { - printf("Problem with forecast resultfile: %s\n", filerespop); - fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop); - } - printf("Computing forecasting: result on file '%s' \n", filerespop); - fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop); +/* strcpy(filerespop,"POP_"); */ +/* strcat(filerespop,fileresu); */ +/* if((ficrespop=fopen(filerespop,"w"))==NULL) { */ +/* printf("Problem with forecast resultfile: %s\n", filerespop); */ +/* fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop); */ +/* } */ +/* printf("Computing forecasting: result on file '%s' \n", filerespop); */ +/* fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop); */ - if (cptcoveff==0) ncodemax[cptcoveff]=1; +/* if (cptcoveff==0) ncodemax[cptcoveff]=1; */ - /* if (mobilav!=0) { */ - /* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */ - /* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ */ - /* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */ - /* printf(" Error in movingaverage mobilav=%d\n",mobilav); */ - /* } */ - /* } */ +/* /\* if (mobilav!=0) { *\/ */ +/* /\* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */ +/* /\* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ *\/ */ +/* /\* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); *\/ */ +/* /\* printf(" Error in movingaverage mobilav=%d\n",mobilav); *\/ */ +/* /\* } *\/ */ +/* /\* } *\/ */ - stepsize=(int) (stepm+YEARM-1)/YEARM; - if (stepm<=12) stepsize=1; +/* stepsize=(int) (stepm+YEARM-1)/YEARM; */ +/* if (stepm<=12) stepsize=1; */ - agelim=AGESUP; +/* agelim=AGESUP; */ - hstepm=1; - hstepm=hstepm/stepm; +/* hstepm=1; */ +/* hstepm=hstepm/stepm; */ - if (popforecast==1) { - if((ficpop=fopen(popfile,"r"))==NULL) { - printf("Problem with population file : %s\n",popfile);exit(0); - fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0); - } - popage=ivector(0,AGESUP); - popeffectif=vector(0,AGESUP); - popcount=vector(0,AGESUP); +/* if (popforecast==1) { */ +/* if((ficpop=fopen(popfile,"r"))==NULL) { */ +/* printf("Problem with population file : %s\n",popfile);exit(0); */ +/* fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0); */ +/* } */ +/* popage=ivector(0,AGESUP); */ +/* popeffectif=vector(0,AGESUP); */ +/* popcount=vector(0,AGESUP); */ - i=1; - while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1; +/* i=1; */ +/* while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1; */ - imx=i; - for (i=1; i=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ - nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); - nhstepm = nhstepm/hstepm; +/* for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ */ +/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); */ +/* nhstepm = nhstepm/hstepm; */ - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - oldm=oldms;savm=savms; - hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); +/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */ +/* oldm=oldms;savm=savms; */ +/* hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */ - for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedatem+YEARM*cpt)) { - fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); - } - for(j=1; j<=nlstate+ndeath;j++) { - kk1=0.;kk2=0; - for(i=1; i<=nlstate;i++) { - if (mobilav==1) - kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; - else { - kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; - } - } - if (h==(int)(calagedatem+12*cpt)){ - tabpop[(int)(agedeb)][j][cptcod]=kk1; - /*fprintf(ficrespop," %.3f", kk1); - if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/ - } - } - for(i=1; i<=nlstate;i++){ - kk1=0.; - for(j=1; j<=nlstate;j++){ - kk1= kk1+tabpop[(int)(agedeb)][j][cptcod]; - } - tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)]; - } +/* for (h=0; h<=nhstepm; h++){ */ +/* if (h==(int) (calagedatem+YEARM*cpt)) { */ +/* fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); */ +/* } */ +/* for(j=1; j<=nlstate+ndeath;j++) { */ +/* kk1=0.;kk2=0; */ +/* for(i=1; i<=nlstate;i++) { */ +/* if (mobilav==1) */ +/* kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; */ +/* else { */ +/* kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; */ +/* } */ +/* } */ +/* if (h==(int)(calagedatem+12*cpt)){ */ +/* tabpop[(int)(agedeb)][j][cptcod]=kk1; */ +/* /\*fprintf(ficrespop," %.3f", kk1); */ +/* if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*\/ */ +/* } */ +/* } */ +/* for(i=1; i<=nlstate;i++){ */ +/* kk1=0.; */ +/* for(j=1; j<=nlstate;j++){ */ +/* kk1= kk1+tabpop[(int)(agedeb)][j][cptcod]; */ +/* } */ +/* tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)]; */ +/* } */ - if (h==(int)(calagedatem+12*cpt)) - for(j=1; j<=nlstate;j++) - fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]); - } - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - } - } +/* if (h==(int)(calagedatem+12*cpt)) */ +/* for(j=1; j<=nlstate;j++) */ +/* fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]); */ +/* } */ +/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */ +/* } */ +/* } */ - /******/ +/* /\******\/ */ - for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) { - fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); - for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ - nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); - nhstepm = nhstepm/hstepm; +/* for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) { */ +/* fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); */ +/* for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ */ +/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); */ +/* nhstepm = nhstepm/hstepm; */ - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - oldm=oldms;savm=savms; - hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); - for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedatem+YEARM*cpt)) { - fprintf(ficresf,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); - } - for(j=1; j<=nlstate+ndeath;j++) { - kk1=0.;kk2=0; - for(i=1; i<=nlstate;i++) { - kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod]; - } - if (h==(int)(calagedatem+12*cpt)) fprintf(ficresf," %15.2f", kk1); - } - } - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - } - } - } - } +/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */ +/* oldm=oldms;savm=savms; */ +/* hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */ +/* for (h=0; h<=nhstepm; h++){ */ +/* if (h==(int) (calagedatem+YEARM*cpt)) { */ +/* fprintf(ficresf,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); */ +/* } */ +/* for(j=1; j<=nlstate+ndeath;j++) { */ +/* kk1=0.;kk2=0; */ +/* for(i=1; i<=nlstate;i++) { */ +/* kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod]; */ +/* } */ +/* if (h==(int)(calagedatem+12*cpt)) fprintf(ficresf," %15.2f", kk1); */ +/* } */ +/* } */ +/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */ +/* } */ +/* } */ +/* } */ +/* } */ - /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */ +/* /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */ - if (popforecast==1) { - free_ivector(popage,0,AGESUP); - free_vector(popeffectif,0,AGESUP); - free_vector(popcount,0,AGESUP); - } - free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - fclose(ficrespop); -} /* End of popforecast */ +/* if (popforecast==1) { */ +/* free_ivector(popage,0,AGESUP); */ +/* free_vector(popeffectif,0,AGESUP); */ +/* free_vector(popcount,0,AGESUP); */ +/* } */ +/* free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */ +/* free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */ +/* fclose(ficrespop); */ +/* } /\* End of popforecast *\/ */ int fileappend(FILE *fichier, char *optionfich) { @@ -7444,85 +7609,87 @@ int readdata(char datafile[], int firsto /* Loops on waves */ for (j=maxwav;j>=1;j--){ for (iv=nqtv;iv>=1;iv--){ /* Loop on time varying quantitative variables */ - cutv(stra, strb, line, ' '); - if(strb[0]=='.') { /* Missing value */ - lval=-1; - cotqvar[j][iv][i]=-1; /* 0.0/0.0 */ - if(isalpha(strb[1])) { /* .m or .d Really Missing value */ - printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j); - fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog); - return 1; - } - }else{ - errno=0; - /* what_kind_of_number(strb); */ - dval=strtod(strb,&endptr); - /* if( strb[0]=='\0' || (*endptr != '\0')){ */ - /* if(strb != endptr && *endptr == '\0') */ - /* dval=dlval; */ - /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */ - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav); - fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog); - return 1; - } - cotqvar[j][iv][i]=dval; - } - strcpy(line,stra); + cutv(stra, strb, line, ' '); + if(strb[0]=='.') { /* Missing value */ + lval=-1; + cotqvar[j][iv][i]=-1; /* 0.0/0.0 */ + cotvar[j][ntv+iv][i]=-1; /* For performance reasons */ + if(isalpha(strb[1])) { /* .m or .d Really Missing value */ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog); + return 1; + } + }else{ + errno=0; + /* what_kind_of_number(strb); */ + dval=strtod(strb,&endptr); + /* if( strb[0]=='\0' || (*endptr != '\0')){ */ + /* if(strb != endptr && *endptr == '\0') */ + /* dval=dlval; */ + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog); + return 1; + } + cotqvar[j][iv][i]=dval; + cotvar[j][ntv+iv][i]=dval; + } + strcpy(line,stra); }/* end loop ntqv */ for (iv=ntv;iv>=1;iv--){ /* Loop on time varying dummies */ - cutv(stra, strb, line, ' '); - if(strb[0]=='.') { /* Missing value */ - lval=-1; - }else{ - errno=0; - lval=strtol(strb,&endptr,10); - /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav); - fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog); - return 1; - } - } - if(lval <-1 || lval >1){ - printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + cutv(stra, strb, line, ' '); + if(strb[0]=='.') { /* Missing value */ + lval=-1; + }else{ + errno=0; + lval=strtol(strb,&endptr,10); + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog); + return 1; + } + } + if(lval <-1 || lval >1){ + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ - For example, for multinomial values like 1, 2 and 3,\n \ - build V1=0 V2=0 for the reference value (1),\n \ - V1=1 V2=0 for (2) \n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ - output of IMaCh is often meaningless.\n \ + output of IMaCh is often meaningless.\n \ Exiting.\n",lval,linei, i,line,j); - fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ - For example, for multinomial values like 1, 2 and 3,\n \ - build V1=0 V2=0 for the reference value (1),\n \ - V1=1 V2=0 for (2) \n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ - output of IMaCh is often meaningless.\n \ + output of IMaCh is often meaningless.\n \ Exiting.\n",lval,linei, i,line,j);fflush(ficlog); - return 1; - } - cotvar[j][iv][i]=(double)(lval); - strcpy(line,stra); + return 1; + } + cotvar[j][iv][i]=(double)(lval); + strcpy(line,stra); }/* end loop ntv */ /* Statuses at wave */ cutv(stra, strb, line, ' '); if(strb[0]=='.') { /* Missing value */ - lval=-1; + lval=-1; }else{ - errno=0; - lval=strtol(strb,&endptr,10); - /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); - fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog); - return 1; - } + errno=0; + lval=strtol(strb,&endptr,10); + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog); + return 1; + } } s[j][i]=lval; @@ -7680,16 +7847,58 @@ int readdata(char datafile[], int firsto return (1); } -void removespace(char *str) { - char *p1 = str, *p2 = str; +void removespace(char **stri){/*, char stro[]) {*/ + char *p1 = *stri, *p2 = *stri; do while (*p2 == ' ') p2++; while (*p1++ == *p2++); + *stri=p1; } -int decodemodel ( char model[], int lastobs) - /**< This routine decode the model and returns: +int decoderesult ( char resultline[]) +/**< This routine decode one result line and returns the combination # of dummy covariates only **/ +{ + int j=0, k=0; + char resultsav[MAXLINE]; + char stra[80], strb[80], strc[80], strd[80],stre[80]; + + removespace(&resultline); + printf("decoderesult:%s\n",resultline); + + if (strstr(resultline,"v") !=0){ + printf("Error. 'v' must be in upper case 'V' result: %s ",resultline); + fprintf(ficlog,"Error. 'v' must be in upper case result: %s ",resultline);fflush(ficlog); + return 1; + } + trimbb(resultsav, resultline); + if (strlen(resultsav) >1){ + j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */ + } + + for(k=1; k<=j;k++){ /* Loop on total covariates of the model */ + 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 */ + cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */ + Tvalsel[k]=atof(strc); /* 1 */ + + cutl(strc,stre,strd,'V'); /* strd='V4' strc=4 stre='V' */; + Tvarsel[k]=atoi(strc); + /* Typevarsel[k]=1; /\* 1 for age product *\/ */ + /* cptcovsel++; */ + if (nbocc(stra,'=') >0) + strcpy(resultsav,stra); /* and analyzes it */ + } + return (0); +} +int selected( int kvar){ /* Selected combination of covariates */ + if(Tvarsel[kvar]) + return (0); + else + return(1); +} +int decodemodel( char model[], int lastobs) + /**< This routine decodes the model and returns: * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age * - nagesqr = 1 if age*age in the model, otherwise 0. * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age @@ -7705,7 +7914,7 @@ int decodemodel ( char model[], int last */ { int i, j, k, ks; - int j1, k1, k2; + int j1, k1, k2, k3, k4; char modelsav[80]; char stra[80], strb[80], strc[80], strd[80],stre[80]; char *strpt; @@ -7727,21 +7936,21 @@ int decodemodel ( char model[], int last 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){ @@ -7806,70 +8015,72 @@ int decodemodel ( char model[], int last /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */ /* * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ - for(k=cptcovt; k>=1;k--) /**< Number of covariates not including constant and age, neither age*age*/ - Tvar[k]=0; + for(k=cptcovt; k>=1;k--){ /**< Number of covariates not including constant and age, neither age*age*/ + Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0; + } cptcovage=0; for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */ - cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' + cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ - if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ - /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ - /*scanf("%d",i);*/ - if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ - cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ - if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ - /* covar is not filled and then is empty */ - cptcovprod--; - cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ - Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ - Typevar[k]=1; /* 1 for age product */ - cptcovage++; /* Sums the number of covariates which include age as a product */ - Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ - /*printf("stre=%s ", stre);*/ - } else if (strcmp(strd,"age")==0) { /* or age*Vn */ - cptcovprod--; - cutl(stre,strb,strc,'V'); - Tvar[k]=atoi(stre); - Typevar[k]=1; /* 1 for age product */ - cptcovage++; - Tage[cptcovage]=k; - } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ - /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ - cptcovn++; - cptcovprodnoage++;k1++; - cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ - Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but - because this model-covariate is a construction we invent a new column - which is after existing variables ncovcol+nqv+ntv+nqtv + k1 - If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 - Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ - Typevar[k]=2; /* 2 for double fixed dummy covariates */ - cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ - Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ - Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ - Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ - k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */ - /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */ - /* Tvar[cptcovt+k2+1]=Tvard[k1][2]; /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */ + if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ + /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ + /*scanf("%d",i);*/ + if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ + cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ + if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ + /* covar is not filled and then is empty */ + cptcovprod--; + cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ + Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ + Typevar[k]=1; /* 1 for age product */ + cptcovage++; /* Sums the number of covariates which include age as a product */ + Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ + /*printf("stre=%s ", stre);*/ + } else if (strcmp(strd,"age")==0) { /* or age*Vn */ + cptcovprod--; + cutl(stre,strb,strc,'V'); + Tvar[k]=atoi(stre); + Typevar[k]=1; /* 1 for age product */ + cptcovage++; + Tage[cptcovage]=k; + } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ + /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ + cptcovn++; + cptcovprodnoage++;k1++; + cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ + Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but + because this model-covariate is a construction we invent a new column + which is after existing variables ncovcol+nqv+ntv+nqtv + k1 + If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 + Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ + Typevar[k]=2; /* 2 for double fixed dummy covariates */ + cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ + Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ + Tposprod[k]=k1; /* Tpsprod[3]=1, Tposprod[2]=5 */ + Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ + Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ + k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */ + /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */ + /* Tvar[cptcovt+k2+1]=Tvard[k1][2]; /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */ /*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */ - /* 1 2 3 4 5 | Tvar[5+1)=1, Tvar[7]=2 */ - for (i=1; i<=lastobs;i++){ - /* Computes the new covariate which is a product of - covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ - covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; - } - } /* End age is not in the model */ - } /* End if model includes a product */ - else { /* no more sum */ - /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ - /* scanf("%d",i);*/ - cutl(strd,strc,strb,'V'); - ks++; /**< Number of simple covariates*/ - cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */ - Tvar[k]=atoi(strd); - Typevar[k]=0; /* 0 for simple covariates */ - } - strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ + /* 1 2 3 4 5 | Tvar[5+1)=1, Tvar[7]=2 */ + for (i=1; i<=lastobs;i++){ + /* Computes the new covariate which is a product of + covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ + covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; + } + } /* End age is not in the model */ + } /* End if model includes a product */ + else { /* no more sum */ + /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ + /* scanf("%d",i);*/ + cutl(strd,strc,strb,'V'); + ks++; /**< Number of simple covariates dummy or quantitative, fixe or varying */ + cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */ + Tvar[k]=atoi(strd); + Typevar[k]=0; /* 0 for simple covariates */ + } + strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); scanf("%d",i);*/ } /* end of loop + on total covariates */ @@ -7885,116 +8096,280 @@ int decodemodel ( char model[], int last scanf("%d ",i);*/ -/* Decodemodel knows only the grammar (simple, product, age*) of the model but not what kind - of variable (dummy vs quantitative, fixed vs time varying) is behind */ +/* Until here, decodemodel knows only the grammar (simple, product, age*) of the model but not what kind + of variable (dummy vs quantitative, fixed vs time varying) is behind. But we know the # of each. */ /* ncovcol= 1, nqv=1 | ntv=2, nqtv= 1 = 5 possible variables data: 2 fixed 3, varying model= V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place k = 1 2 3 4 5 6 7 8 9 Tvar[k]= 5 4 3 1+1+2+1+1=6 5 2 7 1 5 Typevar[k]= 0 0 0 2 1 0 2 1 1 - Fixed[Tvar[k]]1 1 1 1 2 0 1 2 3 - Dummy[Tvar[k]]1 0 0 0 2 1 1 2 3 + Fixed[k] 1 1 1 1 3 0 0 or 2 2 3 + Dummy[k] 1 0 0 0 3 1 1 2 3 + Tmodelind[combination of covar]=k; */ /* Dispatching between quantitative and time varying covariates */ /* If Tvar[k] >ncovcol it is a product */ /* Tvar[k] is the value n of Vn with n varying for 1 to nvcol, or p Vp=Vn*Vm for product */ /* Computing effective variables, ie used by the model, that is from the cptcovt variables */ - 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 */ - Fixed[Tvar[k]]= 0; - Dummy[Tvar[k]]= 0; + printf("Model=%s\n\ +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); + fprintf(ficlog,"Model=%s\n\ +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, ncovf=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 || Typevar[k]==2)){ /* Simple or product fixed dummy (<=ncovcol) covariates */ + Fixed[k]= 0; + Dummy[k]= 0; ncoveff++; - }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ - Fixed[Tvar[k]]= 0; - Dummy[Tvar[k]]= 1; - nqfveff++; /* Only simple fixed quantitative variable */ + ncovf++; + modell[k].maintype= FTYPE; + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + 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++; + modell[k].maintype= FTYPE; + modell[k].subtype= FQ; + ncovf++; + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + 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[Tvar[k]]= 1; - Dummy[Tvar[k]]= 0; + Fixed[k]= 1; + Dummy[k]= 0; ntveff++; /* Only simple time varying dummy variable */ - }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){ - if( Typevar[k]==0){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - nqtveff++;/* Only simple time varying quantitative variable */ - } - }else if (Typevar[k] == 2) { - for(k1=1; k1 <= cptcovprodnoage; k1++){ - if(Tvard[k1][1] <=ncovcol){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 0; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[Tvar[k]]= 0; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 0; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - } - }else if(Tvard[k1][1] <=ncovcol+nqv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[Tvar[k]]= 0; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[Tvar[k]]= 0; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - } - }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 0; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - } - }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ - if(Tvard[k1][2] <=ncovcol){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ - Fixed[Tvar[k]]= 1; - Dummy[Tvar[k]]= 1; - } - }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]); - } + modell[k].maintype= VTYPE; + modell[k].subtype= VD; + ncovv++; /* Only simple time varying variables */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + 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){ /* Only simple time varying quantitative variable V5*/ + Fixed[k]= 1; + Dummy[k]= 1; + nqtveff++; + modell[k].maintype= VTYPE; + modell[k].subtype= VQ; + ncovv++; /* Only simple time varying variables */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + 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 */ + ncova++; + TvarA[ncova]=Tvar[k]; + TvarAind[ncova]=k; + 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; + 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; + 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; + 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; + modell[k].maintype= FTYPE; + modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */ + ncovf++; /* Fixed variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }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 */ + ncovf++; /* Varying variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }else if(Tvard[k1][1] <=ncovcol+nqv){ + 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 */ + ncovf++; /* Fixed variables without age */ + TvarF[ncovf]=Tvar[k]; + TvarFind[ncovf]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ + 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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ + 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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + }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 */ + ncovv++; /* Varying variables without age */ + TvarV[ncovv]=Tvar[k]; + TvarVind[ncovv]=k; + } + }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]); } /* 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[Tvar[k]],Dummy[Tvar[k]]); - fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[Tvar[k]],Dummy[Tvar[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 0){ /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */ /* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */ - bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k); + bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k); }else if (mobilavproj == 0){ - printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj); - fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj); - exit(1); + printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj); + fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj); + exit(1); }else{ - /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */ - bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k); + /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */ + bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k); } fprintf(ficresplb,"%.0f ",age ); - for(j=1;j<=nqfveff;j++) - fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + for(j=1;j<=cptcoveff;j++) + fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); tot=0.; for(i=1; i<=nlstate;i++){ - tot += bprlim[i][i]; - fprintf(ficresplb," %.5f", bprlim[i][i]); + tot += bprlim[i][i]; + fprintf(ficresplb," %.5f", bprlim[i][i]); } fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp); } /* Age */ @@ -8519,13 +8890,13 @@ int hPijx(double *p, int bage, int fage) /* hstepm=1; aff par mois*/ pstamp(ficrespij); fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); - i1= pow(2,nqfveff); + i1= pow(2,cptcoveff); /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */ /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */ /* k=k+1; */ - for (k=1; k <= (int) pow(2,nqfveff); k++){ + for (k=1; k <= (int) pow(2,cptcoveff); k++){ fprintf(ficrespij,"\n#****** "); - for(j=1;j<=nqfveff;j++) + for(j=1;j<=cptcoveff;j++) fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficrespij,"******\n"); @@ -8591,13 +8962,13 @@ int hPijx(double *p, int bage, int fage) /* hstepm=1; aff par mois*/ pstamp(ficrespijb); fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x "); - i1= pow(2,nqfveff); + i1= pow(2,cptcoveff); /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */ /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */ /* k=k+1; */ - for (k=1; k <= (int) pow(2,nqfveff); k++){ + for (k=1; k <= (int) pow(2,cptcoveff); k++){ fprintf(ficrespijb,"\n#****** "); - for(j=1;j<=nqfveff;j++) + for(j=1;j<=cptcoveff;j++) fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficrespijb,"******\n"); if(invalidvarcomb[k]){ @@ -8682,6 +9053,8 @@ int main(int argc, char *argv[]) char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; char model[MAXLINE], modeltemp[MAXLINE]; + char resultline[MAXLINE]; + char pathr[MAXLINE], pathimach[MAXLINE]; char *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; @@ -8993,7 +9366,7 @@ int main(int argc, char *argv[]) covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */ - cotvar=ma3x(1,maxwav,1,ntv,1,n); /**< Time varying covariate */ + cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/ cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */ 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 @@ -9236,7 +9609,25 @@ Please run with mle=-1 to get a correct k=2 V1 Tvar[k=2]= 1 (from V1) k=1 Tvar[1]=2 (from V2) */ - Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ + + Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ + TvarF=ivector(1,NCOVMAX); /* */ + TvarFind=ivector(1,NCOVMAX); /* */ + TvarV=ivector(1,NCOVMAX); /* */ + TvarVind=ivector(1,NCOVMAX); /* */ + TvarA=ivector(1,NCOVMAX); /* */ + TvarAind=ivector(1,NCOVMAX); /* */ + TvarFD=ivector(1,NCOVMAX); /* */ + TvarFDind=ivector(1,NCOVMAX); /* */ + TvarFQ=ivector(1,NCOVMAX); /* */ + TvarFQind=ivector(1,NCOVMAX); /* */ + TvarVD=ivector(1,NCOVMAX); /* */ + TvarVDind=ivector(1,NCOVMAX); /* */ + TvarVQ=ivector(1,NCOVMAX); /* */ + TvarVQind=ivector(1,NCOVMAX); /* */ + + Tvalsel=vector(1,NCOVMAX); /* */ + Tvarsel=ivector(1,NCOVMAX); /* */ Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */ Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */ Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */ @@ -9249,9 +9640,11 @@ Please run with mle=-1 to get a correct ncovcol + k1 If already ncovcol=4 and model=V2+V1+V1*V4+age*V3 Tvar[3=V1*V4]=4+1 etc */ - Tprod=ivector(1,NCOVMAX); /* Gives the position of a product */ + Tprod=ivector(1,NCOVMAX); /* Gives the k position of the k1 product */ + Tposprod=ivector(1,NCOVMAX); /* Gives the k1 product from the k position */ /* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3 if V2+V1+V1*V4+age*V3+V3*V2 TProd[k1=2]=5 (V3*V2) + Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2 */ Tvaraff=ivector(1,NCOVMAX); /* Unclear */ Tvard=imatrix(1,NCOVMAX,1,2); /* n=Tvard[k1][1] and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm @@ -9261,7 +9654,19 @@ Please run with mle=-1 to get a correct 4 covariates (3 plus signs) Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 */ - + Tmodelind=ivector(1,NCOVMAX);/** gives the k model position of an + * individual dummy, fixed or varying: + * Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4, + * 3, 1, 0, 0, 0, 0, 0, 0}, + * model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 , + * V1 df, V2 qf, V3 & V4 dv, V5 qv + * Tmodelind[1]@9={9,0,3,2,}*/ + TmodelInvind=ivector(1,NCOVMAX); /* TmodelInvind=Tvar[k]- ncovcol-nqv={5-2-1=2,*/ + TmodelInvQind=ivector(1,NCOVMAX);/** gives the k model position of an + * individual quantitative, fixed or varying: + * Tmodelqind[1]=1,Tvaraff[1]@9={4, + * 3, 1, 0, 0, 0, 0, 0, 0}, + * model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ /* Main decodemodel */ @@ -9323,17 +9728,17 @@ Please run with mle=-1 to get a correct cptcoveff=0; if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */ tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ - } - - ncovcombmax=pow(2,cptcoveff); - invalidvarcomb=ivector(1, ncovcombmax); - for(i=1;iDatafile=%s Firstpass=%d La /* Calculates basic frequencies. Computes observed prevalence at single age and for any valid combination of covariates and prints on file fileres'p'. */ - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ - firstpass, lastpass, stepm, weightopt, model); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ + firstpass, lastpass, stepm, weightopt, model); fprintf(fichtm,"\n"); fprintf(fichtm,"
Total number of observations=%d
\n\ @@ -10027,16 +10432,54 @@ Please run with mle=-1 to get a correct fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); /* day and month of proj2 are not used but only year anproj2.*/ + /* Results */ + while(fgets(line, MAXLINE, ficpar)) { + /* If line starts with a # it is a comment */ + if (line[0] == '#') { + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + continue; + }else + break; + } + while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){ + if (num_filled == 0) + resultline[0]='\0'; + else if (num_filled != 1){ + printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line); + } + printf("Result %d: result line should be at minimum 'line=' %s, result=%s\n",num_filled, line, resultline); + decoderesult(resultline); + while(fgets(line, MAXLINE, ficpar)) { + /* If line starts with a # it is a comment */ + if (line[0] == '#') { + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + continue; + }else + break; + } + if (feof(ficpar)) + break; + else{ /* Processess output results for this combination of covariate values */ + } + } + + - /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */ + /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */ /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */ replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){ - printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); - fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); }else{ @@ -10074,12 +10517,12 @@ Please run with mle=-1 to get a correct /*#include "hpijx.h"*/ hPijx(p, bage, fage); fclose(ficrespij); - + /* ncovcombmax= pow(2,cptcoveff); */ /*-------------- Variance of one-step probabilities---*/ k=1; varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart); - + /* Prevalence for each covariates in probs[age][status][cov] */ probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); for(i=1;i<=AGESUP;i++) @@ -10089,27 +10532,27 @@ Please run with mle=-1 to get a correct prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); if (mobilav!=0 ||mobilavproj !=0 ) { mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); - for(i=1;i<=AGESUP;i++) - for(j=1;j<=nlstate;j++) - for(k=1;k<=ncovcombmax;k++) - mobaverages[i][j][k]=0.; + for(i=1;i<=AGESUP;i++) + for(j=1;j<=nlstate;j++) + for(k=1;k<=ncovcombmax;k++) + mobaverages[i][j][k]=0.; mobaverage=mobaverages; if (mobilav!=0) { - if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){ - fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); - printf(" Error in movingaverage mobilav=%d\n",mobilav); - } + if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } /* /\* Prevalence for each covariates in probs[age][status][cov] *\/ */ /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */ else if (mobilavproj !=0) { - if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){ - fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj); - printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj); - } + if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){ + fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj); + printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj); + } } }/* end if moving average */ - + /*---------- Forecasting ------------------*/ /*if((stepm == 1) && (strcmp(model,".")==0)){*/ if(prevfcast==1){ @@ -10158,10 +10601,10 @@ Please run with mle=-1 to get a correct printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout); fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog); - for (k=1; k <= (int) pow(2,cptcoveff); k++){ + for (k=1; k <= (int) pow(2,cptcoveff); k++){ /* For any combination of dummy covariates, fixed and varying */ fprintf(ficreseij,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } fprintf(ficreseij,"******\n"); @@ -10175,7 +10618,7 @@ Please run with mle=-1 to get a correct printf("done evsij\n");fflush(stdout); fprintf(ficlog,"done evsij\n");fflush(ficlog); - /*---------- Health expectancies and variances ------------*/ + /*---------- State-specific expectancies and variances ------------*/ strcpy(filerest,"T_"); @@ -10191,20 +10634,20 @@ Please run with mle=-1 to get a correct strcpy(fileresstde,"STDE_"); strcat(fileresstde,fileresu); if((ficresstdeij=fopen(fileresstde,"w"))==NULL) { - printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); - fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); + printf("Problem with State specific Exp. and std errors resultfile: %s\n", fileresstde); exit(0); + fprintf(ficlog,"Problem with State specific Exp. and std errors resultfile: %s\n", fileresstde); exit(0); } - printf(" Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); - fprintf(ficlog," Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); + printf(" Computing State-specific Expectancies and standard errors: result on file '%s' \n", fileresstde); + fprintf(ficlog," Computing State-specific Expectancies and standard errors: result on file '%s' \n", fileresstde); strcpy(filerescve,"CVE_"); strcat(filerescve,fileresu); if((ficrescveij=fopen(filerescve,"w"))==NULL) { - printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); - fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); + printf("Problem with Covar. State-specific Exp. resultfile: %s\n", filerescve); exit(0); + fprintf(ficlog,"Problem with Covar. State-specific Exp. resultfile: %s\n", filerescve); exit(0); } - printf(" Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); - fprintf(ficlog," Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); + printf(" Computing Covar. of State-specific Expectancies: result on file '%s' \n", filerescve); + fprintf(ficlog," Computing Covar. of State-specific Expectancies: result on file '%s' \n", filerescve); strcpy(fileresv,"V_"); strcat(fileresv,fileresu); @@ -10212,36 +10655,43 @@ Please run with mle=-1 to get a correct printf("Problem with variance resultfile: %s\n", fileresv);exit(0); fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); } - printf(" Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(stdout); - fprintf(ficlog," Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(ficlog); + printf(" Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout); + fprintf(ficlog," Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog); /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ for (k=1; k <= (int) pow(2,cptcoveff); k++){ + printf("\n#****** "); fprintf(ficrest,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficlog,"\n#****** "); + for(j=1;j<=cptcoveff;j++){ + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + } fprintf(ficrest,"******\n"); + fprintf(ficlog,"******\n"); + printf("******\n"); fprintf(ficresstdeij,"\n#****** "); fprintf(ficrescveij,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } fprintf(ficresstdeij,"******\n"); fprintf(ficrescveij,"******\n"); fprintf(ficresvij,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficresvij,"******\n"); eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - printf(" cvevsij %d, ",k); - fprintf(ficlog, " cvevsij %d, ",k); + printf(" cvevsij combination#=%d, ",k); + fprintf(ficlog, " cvevsij combination#=%d, ",k); cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); printf(" end cvevsij \n "); fprintf(ficlog, " end cvevsij \n "); @@ -10255,57 +10705,57 @@ Please run with mle=-1 to get a correct for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ - oldm=oldms;savm=savms; /* ZZ Segmentation fault */ - cptcod= 0; /* To be deleted */ - printf("varevsij %d \n",vpopbased); - fprintf(ficlog, "varevsij %d \n",vpopbased); - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ - fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); - 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,"# 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"); */ - epj=vector(1,nlstate+1); - 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"); - for(age=bage; age <=fage ;age++){ - prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */ - if (vpopbased==1) { - if(mobilav ==0){ - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][k]; - }else{ /* mobilav */ - for(i=1; i<=nlstate;i++) - prlim[i][i]=mobaverage[(int)age][i][k]; - } - } + oldm=oldms;savm=savms; /* ZZ Segmentation fault */ + cptcod= 0; /* To be deleted */ + printf("varevsij vpopbased=%d \n",vpopbased); + fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased); + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ + fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); + 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,"# 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"); */ + epj=vector(1,nlstate+1); + 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"); + for(age=bage; age <=fage ;age++){ + prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */ + if (vpopbased==1) { + if(mobilav ==0){ + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][k]; + }else{ /* mobilav */ + for(i=1; i<=nlstate;i++) + prlim[i][i]=mobaverage[(int)age][i][k]; + } + } - fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); - /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */ - /* printf(" age %4.0f ",age); */ - for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ - for(i=1, epj[j]=0.;i <=nlstate;i++) { - epj[j] += prlim[i][i]*eij[i][j][(int)age]; - /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ - /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ - } - epj[nlstate+1] +=epj[j]; - } - /* printf(" age %4.0f \n",age); */ + fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); + /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */ + /* printf(" age %4.0f ",age); */ + for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ + for(i=1, epj[j]=0.;i <=nlstate;i++) { + epj[j] += prlim[i][i]*eij[i][j][(int)age]; + /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ + /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ + } + epj[nlstate+1] +=epj[j]; + } + /* printf(" age %4.0f \n",age); */ - for(i=1, vepp=0.;i <=nlstate;i++) - for(j=1;j <=nlstate;j++) - vepp += vareij[i][j][(int)age]; - fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); - for(j=1;j <=nlstate;j++){ - fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); - } - fprintf(ficrest,"\n"); - } + for(i=1, vepp=0.;i <=nlstate;i++) + for(j=1;j <=nlstate;j++) + vepp += vareij[i][j][(int)age]; + fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); + for(j=1;j <=nlstate;j++){ + fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); + } + fprintf(ficrest,"\n"); + } } /* End vpopbased */ free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); @@ -10315,23 +10765,12 @@ Please run with mle=-1 to get a correct /*}*/ } /* End k */ - free_vector(weight,1,n); - 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_ivector(tab,1,NCOVMAX); - fclose(ficresstdeij); - fclose(ficrescveij); - fclose(ficresvij); - fclose(ficrest); - printf("done Health expectancies\n");fflush(stdout); - fprintf(ficlog,"done Health expectancies\n");fflush(ficlog); - fclose(ficpar); - - /*------- Variance of period (stable) prevalence------*/ + printf("done State-specific expectancies\n");fflush(stdout); + fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); + + /*------- Variance of period (stable) prevalence------*/ + strcpy(fileresvpl,"VPL_"); strcat(fileresvpl,fileresu); if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { @@ -10340,27 +10779,48 @@ Please run with mle=-1 to get a correct } 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); - + /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ - + for (k=1; k <= (int) pow(2,cptcoveff); k++){ - fprintf(ficresvpl,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficresvpl,"******\n"); - - varpl=matrix(1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; - varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart); - free_matrix(varpl,1,nlstate,(int) bage, (int)fage); + fprintf(ficresvpl,"\n#****** "); + printf("\n#****** "); + fprintf(ficlog,"\n#****** "); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + } + fprintf(ficresvpl,"******\n"); + printf("******\n"); + fprintf(ficlog,"******\n"); + + varpl=matrix(1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart); + free_matrix(varpl,1,nlstate,(int) bage, (int)fage); /*}*/ } - + fclose(ficresvpl); printf("done variance-covariance of period prevalence\n");fflush(stdout); fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); - + + free_vector(weight,1,n); + 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_ivector(tab,1,NCOVMAX); + fclose(ficresstdeij); + fclose(ficrescveij); + fclose(ficresvij); + fclose(ficrest); + fclose(ficpar); + + /*---------- End : free ----------------*/ if (mobilav!=0 ||mobilavproj !=0) free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */ @@ -10368,38 +10828,58 @@ Please run with mle=-1 to get a correct free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); } /* mle==-3 arrives here for freeing */ - /* endfree:*/ - free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); - free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); - free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); - free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); - free_ma3x(cotvar,1,maxwav,1,ntv,1,n); - free_matrix(coqvar,1,maxwav,1,n); - free_matrix(covar,0,NCOVMAX,1,n); - free_matrix(matcov,1,npar,1,npar); - free_matrix(hess,1,npar,1,npar); - /*free_vector(delti,1,npar);*/ - free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - free_matrix(agev,1,maxwav,1,imx); - free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - - free_ivector(ncodemax,1,NCOVMAX); - free_ivector(ncodemaxwundef,1,NCOVMAX); - free_ivector(Dummy,-1,NCOVMAX); - free_ivector(Fixed,-1,NCOVMAX); - free_ivector(Typevar,-1,NCOVMAX); - free_ivector(Tvar,1,NCOVMAX); - free_ivector(Tprod,1,NCOVMAX); - free_ivector(Tvaraff,1,NCOVMAX); - free_ivector(invalidvarcomb,1,ncovcombmax); - free_ivector(Tage,1,NCOVMAX); - - free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); - /* free_imatrix(codtab,1,100,1,10); */ + /* endfree:*/ + free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); + free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); + free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); + free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); + free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n); + free_matrix(coqvar,1,maxwav,1,n); + free_matrix(covar,0,NCOVMAX,1,n); + free_matrix(matcov,1,npar,1,npar); + free_matrix(hess,1,npar,1,npar); + /*free_vector(delti,1,npar);*/ + free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); + free_matrix(agev,1,maxwav,1,imx); + free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); + + free_ivector(ncodemax,1,NCOVMAX); + free_ivector(ncodemaxwundef,1,NCOVMAX); + free_ivector(Dummy,-1,NCOVMAX); + free_ivector(Fixed,-1,NCOVMAX); + free_ivector(Typevar,-1,NCOVMAX); + free_ivector(Tvar,1,NCOVMAX); + free_ivector(TvarFD,1,NCOVMAX); + free_ivector(TvarFDind,1,NCOVMAX); + free_ivector(TvarF,1,NCOVMAX); + free_ivector(TvarFind,1,NCOVMAX); + free_ivector(TvarV,1,NCOVMAX); + free_ivector(TvarVind,1,NCOVMAX); + free_ivector(TvarA,1,NCOVMAX); + free_ivector(TvarAind,1,NCOVMAX); + free_ivector(TvarFQ,1,NCOVMAX); + free_ivector(TvarFQind,1,NCOVMAX); + free_ivector(TvarVD,1,NCOVMAX); + free_ivector(TvarVDind,1,NCOVMAX); + free_ivector(TvarVQ,1,NCOVMAX); + free_ivector(TvarVQind,1,NCOVMAX); + free_ivector(Tvarsel,1,NCOVMAX); + free_vector(Tvalsel,1,NCOVMAX); + free_ivector(Tposprod,1,NCOVMAX); + free_ivector(Tprod,1,NCOVMAX); + free_ivector(Tvaraff,1,NCOVMAX); + free_ivector(invalidvarcomb,1,ncovcombmax); + free_ivector(Tage,1,NCOVMAX); + free_ivector(Tmodelind,1,NCOVMAX); + free_ivector(TmodelInvind,1,NCOVMAX); + free_ivector(TmodelInvQind,1,NCOVMAX); + + free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); + /* free_imatrix(codtab,1,100,1,10); */ fflush(fichtm); fflush(ficgp); - + if((nberr >0) || (nbwarn>0)){ printf("End of Imach with %d errors and/or %d warnings. Please look at the log file for details.\n",nberr,nbwarn); fprintf(ficlog,"End of Imach with %d errors and/or warnings %d. Please look at the log file for details.\n",nberr,nbwarn); @@ -10417,7 +10897,7 @@ Please run with mle=-1 to get a correct printf("Local time at start %s\nLocal time at end %s",strstart, strtend); fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend); printf("Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout)); - + printf("Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time)); fprintf(ficlog,"Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout)); fprintf(ficlog,"Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time)); @@ -10430,17 +10910,17 @@ Please run with mle=-1 to get a correct fclose(ficgp); fclose(ficlog); /*------ End -----------*/ - - - printf("Before Current directory %s!\n",pathcd); + + + printf("Before Current directory %s!\n",pathcd); #ifdef WIN32 - if (_chdir(pathcd) != 0) - printf("Can't move to directory %s!\n",path); - if(_getcwd(pathcd,MAXLINE) > 0) + if (_chdir(pathcd) != 0) + printf("Can't move to directory %s!\n",path); + if(_getcwd(pathcd,MAXLINE) > 0) #else - if(chdir(pathcd) != 0) - printf("Can't move to directory %s!\n", path); - if (getcwd(pathcd, MAXLINE) > 0) + if(chdir(pathcd) != 0) + printf("Can't move to directory %s!\n", path); + if (getcwd(pathcd, MAXLINE) > 0) #endif printf("Current directory %s!\n",pathcd); /*strcat(plotcmd,CHARSEPARATOR);*/ @@ -10466,7 +10946,7 @@ Please run with mle=-1 to get a correct sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); - + if((outcmd=system(plotcmd)) != 0){ printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); @@ -10494,7 +10974,7 @@ Please run with mle=-1 to get a correct else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } - end: +end: while (z[0] != 'q') { printf("\nType q for exiting: "); fflush(stdout); scanf("%s",z);