--- imach/src/imach.c 2016/08/23 07:40:50 1.233 +++ imach/src/imach.c 2016/08/23 16:51:20 1.234 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.233 2016/08/23 07:40:50 brouard Exp $ +/* $Id: imach.c,v 1.234 2016/08/23 16:51:20 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.234 2016/08/23 16:51:20 brouard + *** empty log message *** + Revision 1.233 2016/08/23 07:40:50 brouard Summary: not working @@ -899,12 +902,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.233 2016/08/23 07:40:50 brouard Exp $ */ +/* $Id: imach.c,v 1.234 2016/08/23 16:51:20 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.233 $ $Date: 2016/08/23 07:40:50 $"; +char fullversion[]="$Revision: 1.234 $ $Date: 2016/08/23 16:51:20 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -920,7 +923,8 @@ int cptcoveff=0; /* Total number of cova 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 nsd=0; /**< Total number of single dummy variables (output) */ +int nsq=0; /**< Total number of single quantitative variables (output) */ 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 */ @@ -978,6 +982,7 @@ char fileresv[FILENAMELENGTH]; FILE *ficresvpl; char fileresvpl[FILENAMELENGTH]; char title[MAXLINE]; +char model[MAXLINE]; /**< The model line */ char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH]; char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH]; char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; @@ -1077,6 +1082,30 @@ double ***cotvar; /* Time varying covari double ***cotqvar; /* Time varying quantitative covariate itqv */ double idx; int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ +/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +/*k 1 2 3 4 5 6 7 8 9 */ +/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */ +/* Tndvar[k] 1 2 3 4 5 */ +/*TDvar 4 3 6 7 1 */ /* For outputs only; combination of dummies fixed or varying */ +/* Tns[k] 1 2 2 4 5 */ /* Number of single cova */ +/* TvarsD[k] 1 2 3 */ /* Number of single dummy cova */ +/* TvarsDind 2 3 9 */ /* position K of single dummy cova */ +/* TvarsQ[k] 1 2 */ /* Number of single quantitative cova */ +/* TvarsQind 1 6 */ /* position K of single quantitative cova */ +/* Tprod[i]=k 4 7 */ +/* Tage[i]=k 5 8 */ +/* */ +/* Type */ +/* V 1 2 3 4 5 */ +/* F F V V V */ +/* D Q D D Q */ +/* */ +int *TvarsD; +int *TvarsDind; +int *TvarsQ; +int *TvarsQind; + +/* int *TDvar; /\**< TDvar[1]=4, TDvarF[2]=3, TDvar[3]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */ 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 */ @@ -1336,7 +1365,7 @@ int nbocc(char *s, char occ) i=0; lg=strlen(s); for(i=0; i<= lg; i++) { - if (s[i] == occ ) j++; + if (s[i] == occ ) j++; } return j; } @@ -2227,87 +2256,87 @@ void powell(double p[], double **xi, int if (directest < 0.0) { /* Then we use it for new direction */ #endif #ifdef DEBUGLINMIN - printf("Before linmin in direction P%d-P0\n",n); - for (j=1;j<=n;j++) { - printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - if(j % ncovmodel == 0){ - printf("\n"); - fprintf(ficlog,"\n"); - } - } + printf("Before linmin in direction P%d-P0\n",n); + for (j=1;j<=n;j++) { + printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ + printf("\n"); + fprintf(ficlog,"\n"); + } + } #endif #ifdef LINMINORIGINAL - linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ + linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ #else - linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ - flatdir[i]=flat; /* Function is vanishing in that direction i */ + linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ + flatdir[i]=flat; /* Function is vanishing in that direction i */ #endif - + #ifdef DEBUGLINMIN - for (j=1;j<=n;j++) { - printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - if(j % ncovmodel == 0){ - printf("\n"); - fprintf(ficlog,"\n"); - } - } + for (j=1;j<=n;j++) { + printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ + printf("\n"); + fprintf(ficlog,"\n"); + } + } #endif - for (j=1;j<=n;j++) { - xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ - xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ - } + for (j=1;j<=n;j++) { + xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ + xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ + } #ifdef LINMINORIGINAL #else - for (j=1, flatd=0;j<=n;j++) { - if(flatdir[j]>0) - flatd++; - } - if(flatd >0){ - printf("%d flat directions\n",flatd); - fprintf(ficlog,"%d flat directions\n",flatd); - for (j=1;j<=n;j++) { - if(flatdir[j]>0){ - printf("%d ",j); - fprintf(ficlog,"%d ",j); - } - } - printf("\n"); - fprintf(ficlog,"\n"); - } + for (j=1, flatd=0;j<=n;j++) { + if(flatdir[j]>0) + flatd++; + } + if(flatd >0){ + printf("%d flat directions\n",flatd); + fprintf(ficlog,"%d flat directions\n",flatd); + for (j=1;j<=n;j++) { + if(flatdir[j]>0){ + printf("%d ",j); + fprintf(ficlog,"%d ",j); + } + } + printf("\n"); + fprintf(ficlog,"\n"); + } #endif - printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); - fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); - + printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + #ifdef DEBUG - printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); - fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); - for(j=1;j<=n;j++){ - printf(" %lf",xit[j]); - fprintf(ficlog," %lf",xit[j]); - } - printf("\n"); - fprintf(ficlog,"\n"); + printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); + fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); + for(j=1;j<=n;j++){ + printf(" %lf",xit[j]); + fprintf(ficlog," %lf",xit[j]); + } + printf("\n"); + fprintf(ficlog,"\n"); #endif } /* end of t or directest negative */ #ifdef POWELLNOF3INFF1TEST #else - } /* end if (fptt < fp) */ + } /* end if (fptt < fp) */ #endif #ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */ - } /*NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done */ + } /*NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done */ #else #endif - } /* loop iteration */ + } /* loop iteration */ } - + /**** Prevalence limit (stable or period prevalence) ****************/ - + double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij) -{ - /* Computes the prevalence limit in each live state at age x and for covariate ij by left multiplying the unit - matrix by transitions matrix until convergence is reached with precision ftolpl */ + { + /* Computes the prevalence limit in each live state at age x and for covariate combiation ij by left multiplying the unit + matrix by transitions matrix until convergence is reached with precision ftolpl */ /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ /* Wx is row vector: population in state 1, population in state 2, population dead */ /* or prevalence in state 1, prevalence in state 2, 0 */ @@ -2325,7 +2354,7 @@ double **prevalim(double **prlim, int nl /* {0.51571254859325999, 0.4842874514067399, */ /* 0.51326036147820708, 0.48673963852179264} */ /* If we start from prlim again, prlim tends to a constant matrix */ - + int i, ii,j,k; double *min, *max, *meandiff, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ @@ -2355,19 +2384,31 @@ double **prevalim(double **prlim, int nl cov[2]=agefin; if(nagesqr==1) cov[3]= agefin*agefin;; - for (k=1; k<=cptcovn;k++) { - /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ - /* Here comes the value of the covariate 'ij' */ - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; - /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */ + for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ + /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; + printf("prevalim ij=%d k=%d TvarsD[%d]=%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); + } + for (k=1; k<=nsq;k++) { /* For single varying covariates only */ + /* Here comes the value of quantitative after renumbering k with single quantitative covariates */ + /* cov[2+nagesqr+TvarsQind[k]]=qselvar[k]; */ + printf("prevalim ij=%d k=%d TvarsQind[%d]=%d \n",ij,k,k,TvarsQind[k]); } /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */ - for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; - for (k=1; k<=cptcovprod;k++) /* Useless */ - /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ + for (k=1; k<=cptcovage;k++){ + if(Dummy[Tvar[Tage[k]]]){ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + } else{ + ; + /* cov[2+nagesqr+Tage[k]]=qselvar[k]; */ + } + printf("prevalim Age ij=%d k=%d Tage[%d]=%d \n",ij,k,k,Tage[k]); + } + for (k=1; k<=cptcovprod;k++){ /* */ + printf("prevalim Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=%d, Tvard[%d][2]=%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; - + } /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ @@ -2517,10 +2558,10 @@ Earliest age to start was %d-%d=%d, ncvl } for(j=1; j<=nlstate; j++){ for(i=1;i<=nlstate;i++){ - /* bprlim[i][j]= newm[i][j]/(1-sumnew); */ - bprlim[i][j]= newm[i][j]; - max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */ - min[i]=FMIN(min[i],bprlim[i][j]); + /* bprlim[i][j]= newm[i][j]/(1-sumnew); */ + bprlim[i][j]= newm[i][j]; + max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */ + min[i]=FMIN(min[i],bprlim[i][j]); } } @@ -2715,81 +2756,81 @@ double **bpmij(double **ps, double *cov, /*double t34;*/ int i,j, nc, ii, jj; - for(i=1; i<= nlstate; i++){ - for(j=1; ji s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ - } - ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ - } - } - - for(i=1; i<= nlstate; i++){ - s1=0; - for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ - ps[i][i]=1./(s1+1.); - /* Computing other pijs */ - for(j=1; ji s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ + } + ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ + } + } + + for(i=1; i<= nlstate; i++){ + s1=0; + for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ + ps[i][i]=1./(s1+1.); + /* Computing other pijs */ + for(j=1; j 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. - */ + 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 */ @@ -4078,74 +4119,74 @@ Title=%s
Datafile=%s Firstpass=%d La for (iind=1; iind<=imx; iind++) { /* For each individual iind */ bool=1; if(anyvaryingduminmodel==0){ /* If All fixed covariates */ - if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ /* for (z1=1; z1<= nqfveff; z1++) { */ /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ /* } */ - for (z1=1; z1<=cptcoveff; z1++) { - /* if(Tvaraff[z1] ==-20){ */ - /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ - /* }else if(Tvaraff[z1] ==-10){ */ - /* /\* sumnew+=coqvar[z1][iind]; *\/ */ - /* }else */ - if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ - /* Tests if this individual iind responded to j1 (V4=1 V3=0) */ - bool=0; - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", - bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), - j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ - /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ - } /* Onlyf fixed */ - } /* end z1 */ - } /* cptcovn > 0 */ + for (z1=1; z1<=cptcoveff; z1++) { + /* if(Tvaraff[z1] ==-20){ */ + /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ + /* }else if(Tvaraff[z1] ==-10){ */ + /* /\* sumnew+=coqvar[z1][iind]; *\/ */ + /* }else */ + if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ + /* Tests if this individual iind responded to j1 (V4=1 V3=0) */ + bool=0; + /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", + bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), + j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ + /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ + } /* Onlyf fixed */ + } /* end z1 */ + } /* cptcovn > 0 */ } /* end any */ if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */ - /* for(m=firstpass; m<=lastpass; m++){ */ - for(mi=1; mi=firstpass && m <=lastpass){ - k2=anint[m][iind]+(mint[m][iind]/12.); - /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ - if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ - if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ - if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ - prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ - if (m1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { - dateintsum=dateintsum+k2; - k2cpt++; - /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ - } - } /* end bool 2 */ - } /* end m */ + /* for(m=firstpass; m<=lastpass; m++){ */ + for(mi=1; mi=firstpass && m <=lastpass){ + k2=anint[m][iind]+(mint[m][iind]/12.); + /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ + if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ + if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ + if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ + prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ + if (m1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { + dateintsum=dateintsum+k2; + k2cpt++; + /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ + } + } /* end bool 2 */ + } /* end m */ } /* end bool */ } /* end iind = 1 to imx */ /* prop[s][age] is feeded for any initial and valid live state as well as @@ -4160,9 +4201,9 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtm, "\n

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

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

\n"); @@ -4184,8 +4225,8 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtmfr,"Age "); for(jk=-1; jk <=nlstate+ndeath; jk++){ for(m=-1; m <=nlstate+ndeath; m++){ - if(jk!=0 && m!=0) - fprintf(ficresphtmfr,"%d%d ",jk,m); + if(jk!=0 && m!=0) + fprintf(ficresphtmfr,"%d%d ",jk,m); } } fprintf(ficresphtmfr, "\n"); @@ -7847,23 +7888,29 @@ int readdata(char datafile[], int firsto return (1); } -void removespace(char **stri){/*, char stro[]) {*/ +void removefirstspace(char **stri){/*, char stro[]) {*/ char *p1 = *stri, *p2 = *stri; - do - while (*p2 == ' ') - p2++; - while (*p1++ == *p2++); - *stri=p1; + if (*p2 == ' ') + p2++; + /* while ((*p1++ = *p2++) !=0) */ + /* ; */ + /* do */ + /* while (*p2 == ' ') */ + /* p2++; */ + /* while (*p1++ == *p2++); */ + *stri=p2; } int decoderesult ( char resultline[]) /**< This routine decode one result line and returns the combination # of dummy covariates only **/ { - int j=0, k=0; + int j=0, k=0, k1=0, k2=0, match=0; char resultsav[MAXLINE]; + int resultmodel[MAXLINE]; + int modelresult[MAXLINE]; char stra[80], strb[80], strc[80], strd[80],stre[80]; - removespace(&resultline); + removefirstspace(&resultline); printf("decoderesult:%s\n",resultline); if (strstr(resultline,"v") !=0){ @@ -7875,13 +7922,19 @@ int decoderesult ( char 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 */ + if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ + printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); + fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); + } + for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */ + if(nbocc(resultsav,'=') >1){ + cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' + resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */ + cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */ + }else + cutl(strc,strd,resultsav,'='); 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 *\/ */ @@ -7889,6 +7942,42 @@ int decoderesult ( char resultline[]) if (nbocc(stra,'=') >0) strcpy(resultsav,stra); /* and analyzes it */ } + /* Checking if no missing or useless values in comparison of current model needs */ + for(k1=1; k1<= cptcovt ;k1++){ /* model line */ + if(Typevar[k1]==0){ + match=0; + for(k2=1; k2 <=j;k2++){ + if(Tvar[k1]==Tvarsel[k2]) { + modelresult[k2]=k1; + match=1; + break; + } + } + if(match == 0){ + printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + } + } + } + + for(k2=1; k2 <=j;k2++){ /* result line */ + match=0; + for(k1=1; k1<= cptcovt ;k1++){ /* model line */ + if(Typevar[k1]==0){ + if(Tvar[k1]==Tvarsel[k2]) { + resultmodel[k1]=k2; + ++match; + } + } + } + if(match == 0){ + printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + }else if(match > 1){ + printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); + } + } + + /* We need to deduce which combination number is chosen and save quantitative values */ + return (0); } int selected( int kvar){ /* Selected combination of covariates */ @@ -7936,21 +8025,21 @@ int decodemodel( char model[], int lasto if ((strpt=strstr(model,"age*age")) !=0){ printf(" strpt=%s, model=%s\n",strpt, model); if(strpt != model){ - printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ corresponding column of parameters.\n",model); - fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ corresponding column of parameters.\n",model); fflush(ficlog); - return 1; + return 1; } nagesqr=1; if (strstr(model,"+age*age") !=0) - substrchaine(modelsav, model, "+age*age"); + substrchaine(modelsav, model, "+age*age"); else if (strstr(model,"age*age+") !=0) - substrchaine(modelsav, model, "age*age+"); + substrchaine(modelsav, model, "age*age+"); else - substrchaine(modelsav, model, "age*age"); + substrchaine(modelsav, model, "age*age"); }else nagesqr=0; if (strlen(modelsav) >1){ @@ -8020,67 +8109,67 @@ int decodemodel( char model[], int lasto } 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 */ - 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) *\/ */ + 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 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 */ + /* 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 */ @@ -8120,221 +8209,243 @@ Typevar: 0 for simple covariate (dummy, 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 */ + for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ + if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ Fixed[k]= 0; Dummy[k]= 0; ncoveff++; ncovf++; - modell[k].maintype= FTYPE; - TvarF[ncovf]=Tvar[k]; - TvarFind[ncovf]=k; + nsd++; + modell[k].maintype= FTYPE; + TvarsD[nsd]=Tvar[k]; + TvarsDind[nsd]=k; + 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 && Typevar[k]==2){ /* Product of fixed dummy (<=ncovcol) covariates */ + Fixed[k]= 0; + Dummy[k]= 0; + ncoveff++; + 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; + modell[k].maintype= FTYPE; + modell[k].subtype= FQ; + nsq++; + TvarsQ[nsq]=Tvar[k]; + TvarsQind[nsq]=k; ncovf++; - TvarF[ncovf]=Tvar[k]; - TvarFind[ncovf]=k; + 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){ + }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying variables */ Fixed[k]= 1; Dummy[k]= 0; ntveff++; /* Only simple time varying dummy variable */ - modell[k].maintype= VTYPE; - modell[k].subtype= VD; - ncovv++; /* Only simple time varying variables */ - TvarV[ncovv]=Tvar[k]; - TvarVind[ncovv]=k; + modell[k].maintype= VTYPE; + modell[k].subtype= VD; + nsd++; + TvarsD[nsd]=Tvar[k]; + TvarsDind[nsd]=k; + 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; + Fixed[k]= 1; + Dummy[k]= 1; + nqtveff++; + modell[k].maintype= VTYPE; + modell[k].subtype= VQ; + ncovv++; /* Only simple time varying variables */ + nsq++; + TvarsQ[nsq]=Tvar[k]; + TvarsQind[nsq]=k; + 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); + 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; + 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++; */ + 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 *\/ */ + 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 *\/ */ + 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 *\/ */ + 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; - } + 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; - } + 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; - } + 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; - } + 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]); + 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]); @@ -8348,28 +8459,28 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( for(k1=1; k1<= cptcovt;k1++){ for(k2=1; k2