--- imach/src/imach.c 2016/09/02 07:25:01 1.245 +++ imach/src/imach.c 2017/03/08 16:02:28 1.255 @@ -1,6 +1,36 @@ -/* $Id: imach.c,v 1.245 2016/09/02 07:25:01 brouard Exp $ +/* $Id: imach.c,v 1.255 2017/03/08 16:02:28 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.255 2017/03/08 16:02:28 brouard + Summary: IMaCh version 0.99r10 bugs in gnuplot fixed + + Revision 1.254 2017/03/08 07:13:00 brouard + Summary: Fixing data parameter line + + Revision 1.253 2016/12/15 11:59:41 brouard + Summary: 0.99 in progress + + Revision 1.252 2016/09/15 21:15:37 brouard + *** empty log message *** + + Revision 1.251 2016/09/15 15:01:13 brouard + Summary: not working + + Revision 1.250 2016/09/08 16:07:27 brouard + Summary: continue + + Revision 1.249 2016/09/07 17:14:18 brouard + Summary: Starting values from frequencies + + Revision 1.248 2016/09/07 14:10:18 brouard + *** empty log message *** + + Revision 1.247 2016/09/02 11:11:21 brouard + *** empty log message *** + + Revision 1.246 2016/09/02 08:49:22 brouard + *** empty log message *** + Revision 1.245 2016/09/02 07:25:01 brouard *** empty log message *** @@ -110,9 +140,7 @@ Author: Nicolas Brouard Revision 1.210 2015/11/18 17:41:20 brouard - Summary: Start working on projected prevalences - - Revision 1.209 2015/11/17 22:12:03 brouard + Summary: Start working on projected prevalences Revision 1.209 2015/11/17 22:12:03 brouard Summary: Adding ftolpl parameter Author: N Brouard @@ -900,7 +928,7 @@ typedef struct { /* #include */ /* #define _(String) gettext (String) */ -#define MAXLINE 1024 /* Was 256. Overflow with 312 with 2 states and 4 covariates. Should be ok */ +#define MAXLINE 2048 /* Was 256 and 1024. Overflow with 312 with 2 states and 4 covariates. Should be ok */ #define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ @@ -937,12 +965,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.245 2016/09/02 07:25:01 brouard Exp $ */ +/* $Id: imach.c,v 1.255 2017/03/08 16:02:28 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.245 $ $Date: 2016/09/02 07:25:01 $"; +char fullversion[]="$Revision: 1.255 $ $Date: 2017/03/08 16:02:28 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -2375,8 +2403,8 @@ void powell(double p[], double **xi, int flatd++; } if(flatd >0){ - printf("%d flat directions\n",flatd); - fprintf(ficlog,"%d flat directions\n",flatd); + printf("%d flat directions: ",flatd); + fprintf(ficlog,"%d flat directions :",flatd); for (j=1;j<=n;j++) { if(flatdir[j]>0){ printf("%d ",j); @@ -2581,6 +2609,7 @@ Earliest age to start was %d-%d=%d, ncvl /* If we start from prlim again, prlim tends to a constant matrix */ int i, ii,j,k; + int first=0; double *min, *max, *meandiff, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ double **out, cov[NCOVMAX+1], **bmij(); @@ -2706,7 +2735,12 @@ Earliest age to start was %d-%d=%d, ncvl } } /* age loop */ /* After some age loop it doesn't converge */ - printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ + if(first){ + first=1; + printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\ +Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + } + fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ free_vector(min,1,nlstate); @@ -3593,10 +3627,11 @@ double funcone( double *x) 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=10 || firstime ==1){ - printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol); - fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol); + printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); + fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); } @@ -4163,17 +4198,84 @@ void pstamp(FILE *fichier) fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); } +int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) { + + /* y=a+bx regression */ + double sumx = 0.0; /* sum of x */ + double sumx2 = 0.0; /* sum of x**2 */ + double sumxy = 0.0; /* sum of x * y */ + double sumy = 0.0; /* sum of y */ + double sumy2 = 0.0; /* sum of y**2 */ + double sume2; /* sum of square or residuals */ + double yhat; + + double denom=0; + int i; + int ne=*no; + + for ( i=ifi, ne=0;i<=ila;i++) { + if(!isfinite(x[i]) || !isfinite(y[i])){ + /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */ + continue; + } + ne=ne+1; + sumx += x[i]; + sumx2 += x[i]*x[i]; + sumxy += x[i] * y[i]; + sumy += y[i]; + sumy2 += y[i]*y[i]; + denom = (ne * sumx2 - sumx*sumx); + /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */ + } + + denom = (ne * sumx2 - sumx*sumx); + if (denom == 0) { + // vertical, slope m is infinity + *b = INFINITY; + *a = 0; + if (r) *r = 0; + return 1; + } + + *b = (ne * sumxy - sumx * sumy) / denom; + *a = (sumy * sumx2 - sumx * sumxy) / denom; + if (r!=NULL) { + *r = (sumxy - sumx * sumy / ne) / /* compute correlation coeff */ + sqrt((sumx2 - sumx*sumx/ne) * + (sumy2 - sumy*sumy/ne)); + } + *no=ne; + for ( i=ifi, ne=0;i<=ila;i++) { + if(!isfinite(x[i]) || !isfinite(y[i])){ + /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */ + continue; + } + ne=ne+1; + yhat = y[i] - *a -*b* x[i]; + sume2 += yhat * yhat ; + + denom = (ne * sumx2 - sumx*sumx); + /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */ + } + *sb = sqrt(sume2/(ne-2)/(sumx2 - sumx * sumx /ne)); + *sa= *sb * sqrt(sumx2/ne); + + return 0; +} + /************ Frequencies ********************/ -void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ +void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \ int firstpass, int lastpass, int stepm, int weightopt, char model[]) -{ /* Some frequencies */ +{ /* Some frequencies as well as proposing some starting values */ - int i, m, jk, j1, bool, z1,j, k, iv; + int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0; int iind=0, iage=0; int mi; /* Effective wave */ int first; double ***freq; /* Frequencies */ + double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */ + int no; double *meanq; double **meanqt; double *pp, **prop, *posprop, *pospropt; @@ -4182,7 +4284,7 @@ void freqsummary(char fileres[], int ia double agebegin, ageend; pp=vector(1,nlstate); - prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); + prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ @@ -4226,14 +4328,15 @@ Title=%s
Datafile=%s Firstpass=%d La } fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate)

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); - freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); + y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); + x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+4+AGEMARGE); j1=0; /* j=ncoveff; /\* Only fixed dummy covariates *\/ */ j=cptcoveff; /* Only dummy covariates of the model */ if (cptcovn<1) {j=1;ncodemax[1]=1;} - first=1; /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: reference=low_education V1=0,V2=0 @@ -4241,277 +4344,413 @@ Title=%s
Datafile=%s Firstpass=%d La high_educ V1=0 V2=1 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 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++) + dateintsum=0; + k2cpt=0; + + if(cptcoveff == 0 ) + nl=1; /* Constant model only */ + else + nl=2; + for (nj = 1; nj <= nl; nj++){ /* nj= 1 constant model, nl number of loops. */ + if(nj==1) + j=0; /* First pass for the constant */ + else + j=cptcoveff; /* Other passes for the covariate values */ + first=1; + 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 (i=1; i<=nlstate; i++) { 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; - 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.; */ - /* } */ - /* } */ - - dateintsum=0; - k2cpt=0; - /* For that combination of covariate j1, we count and print the frequencies in one pass */ - for (iind=1; iind<=imx; iind++) { /* For each individual iind */ - bool=1; - if(anyvaryingduminmodel==0){ /* If All fixed covariates */ - if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - /* 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; mi0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + /* for (z1=1; z1<= nqfveff; z1++) { */ + /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ + /* } */ + for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ + /* if(Tvaraff[z1] ==-20){ */ + /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ + /* }else if(Tvaraff[z1] ==-10){ */ + /* /\* sumnew+=coqvar[z1][iind]; *\/ */ + /* }else */ + if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */ + /* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */ + bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ + /* 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 */ + }/* end j==0 */ + 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) && (j==0)) { + dateintsum=dateintsum+k2; /* on all covariates ?*/ + k2cpt++; + /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ } - } - }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ - /* bool =0 we keep that guy which corresponds to the combination of dummy values */ - if(bool==1){ - /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] - and mw[mi+1][iind]. dh depends on stepm. */ - agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ - ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */ - if(m >=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 (cptcoveff>0){ - fprintf(ficresp, "\n#********** Variable "); - fprintf(ficresphtm, "\n

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

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

\n"); - fprintf(ficresphtmfr, "**********\n"); - fprintf(ficlog, "**********\n"); - } - fprintf(ficresphtm,""); - for(i=1; i<=nlstate;i++) { - fprintf(ficresp, " Age Prev(%d) N(%d) N ",i,i); - fprintf(ficresphtm, "",i,i); - } - fprintf(ficresp, "\n"); - fprintf(ficresphtm, "\n"); - - /* Header of frequency table by age */ - fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); - 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); - } - } - fprintf(ficresphtmfr, "\n"); - - /* For each age */ - for(iage=iagemin; iage <= iagemax+3; iage++){ - fprintf(ficresphtm,""); - if(iage==iagemax+1){ - fprintf(ficlog,"1"); - fprintf(ficresphtmfr," "); - }else if(iage==iagemax+2){ - fprintf(ficlog,"0"); - fprintf(ficresphtmfr," "); - }else if(iage==iagemax+3){ - fprintf(ficlog,"Total"); - fprintf(ficresphtmfr," "); - }else{ - if(first==1){ - first=0; - printf("See log file for details...\n"); + }else{ + bool=1; + }/* 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 (cptcoveff>0 && j!=0){ + printf( "\n#********** Variable "); + fprintf(ficresp, "\n#********** Variable "); + fprintf(ficresphtm, "\n

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

********** Variable "); + fprintf(ficlog, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++){ + if(!FixedV[Tvaraff[z1]]){ + printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficlog, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + }else{ + printf( "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + } + } + printf( "**********\n#"); + fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "**********

\n"); + fprintf(ficresphtmfr, "**********\n"); + fprintf(ficlog, "**********\n"); + } + fprintf(ficresphtm,"
Age%d%d
0
Unknown
Total
"); + for(i=1; i<=nlstate;i++) { + fprintf(ficresp, " Age Prev(%d) N(%d) N ",i,i); + fprintf(ficresphtm, "",i,i); + } + fprintf(ficresp, "\n"); + fprintf(ficresphtm, "\n"); + + /* Header of frequency table by age */ + fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); + 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); } - 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(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][iage]; - if(pp[jk]>=1.e-10){ + fprintf(ficresphtmfr, "\n"); + + /* For each age */ + for(iage=iagemin; iage <= iagemax+3; iage++){ + fprintf(ficresphtm,""); + if(iage==iagemax+1){ + fprintf(ficlog,"1"); + fprintf(ficresphtmfr," "); + }else if(iage==iagemax+2){ + fprintf(ficlog,"0"); + fprintf(ficresphtmfr," "); + }else if(iage==iagemax+3){ + fprintf(ficlog,"Total"); + fprintf(ficresphtmfr," "); + }else{ if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + 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(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); } - 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]; - } /* 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 + 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]; + } /* 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] */ - 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){ + for(jk=1; jk <=nlstate ; jk++){ 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]);*/ + 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); + } } - 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]); + } + /* 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]; + } + fprintf(ficresphtmfr,"\n "); + if(iage <= iagemax){ + fprintf(ficresp,"\n"); + fprintf(ficresphtm,"\n"); } - 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]); + if(first==1) + 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); + }else{ + fprintf(ficresphtm,"",pospropt[jk]/posproptt,pospropt[jk],posproptt); + } + } + fprintf(ficresphtm,"\n"); + fprintf(ficresphtm,"
Age%d%d
%d
0
Unknown
Total
%d%d%.5f%.0f%.0f%d%.5f%.0f%.0f%dNaNq%.0f%.0f%dNaNq%.0f%.0f%.0f
TotNanq%.0f%.0f%.5f%.0f%.0f
\n"); + fprintf(ficresphtmfr,"\n"); + if(posproptt < 1.e-5){ + fprintf(ficresphtm,"\n

This combination (%d) is not valid and no result will be produced

",j1); + fprintf(ficresphtmfr,"\n

This combination (%d) is not valid and no result will be produced

",j1); + fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); + invalidvarcomb[j1]=1; + }else{ + fprintf(ficresphtm,"\n

This combination (%d) is valid and result will be produced.

",j1); + invalidvarcomb[j1]=0; + } + fprintf(ficresphtmfr,"\n"); + fprintf(ficlog,"\n"); + if(j!=0){ + printf("#Freqsummary: Starting values for combination j1=%d:\n", j1); + for(i=1,jk=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */ + if(jj==1){ /* Constant case (in fact cste + age) */ + if(j1==1){ /* All dummy covariates to zero */ + freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */ + freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */ + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3])); + fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); + } + }else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */ + for(iage=iagemin; iage <= iagemax+3; iage++){ + x[iage]= (double)iage; + y[iage]= log(freq[i][k][iage]/freq[i][i][iage]); + /* printf("i=%d, k=%d, jk=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,jk,j1,jj, iage, y[iage]); */ + } + linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */ + pstart[jk]=b; + pstart[jk-1]=a; + }else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ + printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); + printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); + pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])); + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4])); + }else{ /* Other cases, like quantitative fixed or varying covariates */ + ; + } + /* printf("%12.7f )", param[i][jj][k]); */ + /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ + jk++; + } /* end jj */ + } /* end k!= i */ + } /* end k */ + } /* end i, jk */ + } /* end j !=0 */ + } /* end selected combination of covariate j1 */ + if(j==0){ /* We can estimate starting values from the occurences in each case */ + printf("#Freqsummary: Starting values for the constants:\n"); + fprintf(ficlog,"\n"); + for(i=1,jk=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + for(jj=1; jj <=ncovmodel; jj++){ + pstart[jk]=p[jk]; /* Setting pstart to p values by default */ + if(jj==1){ /* Age has to be done */ + pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); + printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + } + /* printf("%12.7f )", param[i][jj][k]); */ + /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ + jk++; } - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); + printf("\n"); + fprintf(ficlog,"\n"); } - if(jk!=0 && m!=0) - fprintf(ficresphtmfr,"%.0f ",freq[jk][m][iage]); } - } /* end loop jk */ - posproptt=0.; - for(jk=1; jk <=nlstate; jk++){ - posproptt += pospropt[jk]; - } - fprintf(ficresphtmfr,"\n "); - if(iage <= iagemax){ - fprintf(ficresp,"\n"); - fprintf(ficresphtm,"\n"); } - if(first==1) - printf("Others in log...\n"); + printf("#Freqsummary\n"); fprintf(ficlog,"\n"); - } /* end loop age iage */ - fprintf(ficresphtm,"Tot"); - for(jk=1; jk <=nlstate ; jk++){ - if(posproptt < 1.e-5){ - fprintf(ficresphtm,"Nanq%.0f%.0f",pospropt[jk],posproptt); - }else{ - fprintf(ficresphtm,"%.5f%.0f%.0f",pospropt[jk]/posproptt,pospropt[jk],posproptt); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + /* param[i]|j][k]= freq[jk][m][iagemax+3] */ + printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); + fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); + /* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */ + /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ + /* fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ + /* } */ + } + } /* end loop jk */ + + printf("\n"); + fprintf(ficlog,"\n"); + } /* end j=0 */ + } /* end j */ + + if(mle == -2){ /* We want to use these values as starting values */ + for(i=1, jk=1; i <=nlstate; i++){ + for(j=1; j <=nlstate+ndeath; j++){ + if(j!=i){ + /*ca[0]= k+'a'-1;ca[1]='\0';*/ + printf("%1d%1d",i,j); + fprintf(ficparo,"%1d%1d",i,j); + for(k=1; k<=ncovmodel;k++){ + /* printf(" %lf",param[i][j][k]); */ + /* fprintf(ficparo," %lf",param[i][j][k]); */ + p[jk]=pstart[jk]; + printf(" %f ",pstart[jk]); + fprintf(ficparo," %f ",pstart[jk]); + jk++; + } + printf("\n"); + fprintf(ficparo,"\n"); + } } } - fprintf(ficresphtm,"\n"); - fprintf(ficresphtm,"\n"); - fprintf(ficresphtmfr,"\n"); - if(posproptt < 1.e-5){ - fprintf(ficresphtm,"\n

This combination (%d) is not valid and no result will be produced

",j1); - fprintf(ficresphtmfr,"\n

This combination (%d) is not valid and no result will be produced

",j1); - fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); - invalidvarcomb[j1]=1; - }else{ - fprintf(ficresphtm,"\n

This combination (%d) is valid and result will be produced.

",j1); - invalidvarcomb[j1]=0; - } - fprintf(ficresphtmfr,"\n"); - } /* end selected combination of covariate j1 */ + } /* end mle=-2 */ dateintmean=dateintsum/k2cpt; fclose(ficresp); @@ -4519,10 +4758,12 @@ Title=%s
Datafile=%s Firstpass=%d La fclose(ficresphtmfr); free_vector(meanq,1,nqfveff); free_matrix(meanqt,1,lastpass,1,nqtveff); - free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE); + free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); + free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); + free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_vector(pospropt,1,nlstate); free_vector(posprop,1,nlstate); - free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE); + free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_vector(pp,1,nlstate); /* End of freqsummary */ } @@ -4549,7 +4790,7 @@ void prevalence(double ***probs, double iagemin= (int) agemin; iagemax= (int) agemax; /*pp=vector(1,nlstate);*/ - prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); + prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ j1=0; @@ -4559,7 +4800,7 @@ void prevalence(double ***probs, double 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++) + for(iage=iagemin-AGEMARGE; iage <= iagemax+4+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)]); */ @@ -4590,7 +4831,7 @@ void prevalence(double ***probs, double 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){ + if((int)agev[m][i] iagemax+4+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); } @@ -4627,7 +4868,7 @@ void prevalence(double ***probs, double /* 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); + free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+4+AGEMARGE); } /* End of prevalence */ /************* Waves Concatenation ***************/ @@ -6170,7 +6411,7 @@ void printinghtml(char fileresu[], char for(nres=1; nres <= nresult; nres++) /* For each resultline */ for(k1=1; k1<=m;k1++){ /* For each combination of covariate */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; /* for(i1=1; i1<=ncodemax[k1];i1++){ */ @@ -6221,13 +6462,13 @@ divided by h: hPij } /* Period (stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\ + fprintf(fichtm,"
\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d some years earlier, knowing that we will be in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\ ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); } if(backcast==1){ /* Period (stable) back prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\ + fprintf(fichtm,"
\n- Convergence to mixed (stable) back prevalence in state %d. Or probability to be in state %d at a younger age, knowing that we will be in state (1 to %d) at different older ages. %s_%d-%d-%d.svg
\ ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); } } @@ -6302,7 +6543,7 @@ See page 'Matrix of variance-covariance for(nres=1; nres <= nresult; nres++){ /* For each resultline */ for(k1=1; k1<=m;k1++){ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; /* for(i1=1; i1<=ncodemax[k1];i1++){ */ jj1++; @@ -6398,10 +6639,10 @@ void printinggnuplot(char fileresu[], ch for (k1=1; k1<= m ; k1 ++){ /* For each valid combination of covariate */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; /* We are interested in selected combination by the resultline */ - printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); + /* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */ fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ @@ -6410,14 +6651,14 @@ void printinggnuplot(char fileresu[], ch /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */ /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */ - printf(" V%d=%d ",Tvaraff[k],vlv); + /* printf(" V%d=%d ",Tvaraff[k],vlv); */ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); } for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ - printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + /* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */ fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); } - printf("\n#\n"); + /* printf("\n#\n"); */ fprintf(ficgp,"\n#\n"); if(invalidvarcomb[k1]){ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); @@ -6480,7 +6721,7 @@ void printinggnuplot(char fileresu[], ch /*2 eme*/ for (k1=1; k1<= m ; k1 ++){ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ @@ -6540,7 +6781,7 @@ void printinggnuplot(char fileresu[], ch /*3eme*/ for (k1=1; k1<= m ; k1 ++){ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; for (cpt=1; cpt<= nlstate ; cpt ++) { @@ -6589,7 +6830,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u /* Survival functions (period) from state i in state j by initial state i */ for (k1=1; k1<=m; k1++){ /* For each covariate and each value */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/ fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt); @@ -6635,7 +6876,7 @@ set ter svg size 640, 480\nunset log y\n /* Survival functions (period) from state i in state j by final state j */ for (k1=1; k1<= m ; k1++){ /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ 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); @@ -6689,9 +6930,9 @@ set ter svg size 640, 480\nunset log y\n /* CV preval stable (period) for each covariate */ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; - for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ + for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */ 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 */ @@ -6715,12 +6956,12 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); k=3; /* Offset */ - for (i=1; i<= nlstate ; i ++){ + for (i=1; i<= nlstate ; i ++){ /* State of origin */ if(i==1) fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); else fprintf(ficgp,", '' "); - l=(nlstate+ndeath)*(i-1)+1; + l=(nlstate+ndeath)*(i-1)+1; /* 1, 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */ 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); @@ -6736,10 +6977,10 @@ set ter svg size 640, 480\nunset log y\n /* CV back preval stable (period) for each covariate */ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; - 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 (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life ending state */ + fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -6761,16 +7002,16 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); k=3; /* Offset */ - for (i=1; i<= nlstate ; i ++){ + for (i=1; i<= nlstate ; i ++){ /* State of origin */ if(i==1) fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_")); else fprintf(ficgp,", '' "); /* l=(nlstate+ndeath)*(i-1)+1; */ - l=(nlstate+ndeath)*(cpt-1)+1; + l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */ /* 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 */ + fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */ /* for (j=2; j<= nlstate ; j ++) */ /* fprintf(ficgp,"+$%d",k+l+j-1); */ /* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */ @@ -6787,7 +7028,7 @@ set ter svg size 640, 480\nunset log y\n for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k1) + if(m != 1 && TKresult[nres]!= k1) continue; 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); @@ -6918,7 +7159,7 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */ for(jk=1; jk <=m; jk++) /* For each combination of covariate */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= jk) + if(m != 1 && TKresult[nres]!= jk) continue; fprintf(ficgp,"# Combination of dummy jk=%d and ",jk); for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ @@ -7280,7 +7521,7 @@ set ter svg size 640, 480\nunset log y\n /* if (h==(int)(YEARM*yearp)){ */ for(nres=1; nres <= nresult; nres++) /* For each resultline */ for(k=1; k<=i1;k++){ - if(TKresult[nres]!= k) + if(i1 != 1 && TKresult[nres]!= k) continue; if(invalidvarcomb[k]){ printf("\nCombination (%d) projection ignored because no cases \n",k); @@ -8231,6 +8472,11 @@ int decoderesult ( char resultline[], in if (strlen(resultsav) >1){ j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */ } + if(j == 0){ /* Resultline but no = */ + TKresult[nres]=0; /* Combination for the nresult and the model */ + return (0); + } + if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); @@ -8830,7 +9076,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( } int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) -{ +{/* Check ages at death */ int i, m; int firstone=0; @@ -9167,7 +9413,7 @@ int prevalence_limit(double *p, double * for(k=1; k<=i1;k++){ /* For each combination k of dummy covariates in the model */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ - if(TKresult[nres]!= k) + if(i1 != 1 && TKresult[nres]!= k) continue; /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */ @@ -9264,7 +9510,7 @@ int back_prevalence_limit(double *p, dou for(nres=1; nres <= nresult; nres++){ /* For each resultline */ for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ - if(TKresult[nres]!= k) + if(i1 != 1 && TKresult[nres]!= k) continue; //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov)); fprintf(ficresplb,"#******"); @@ -9323,6 +9569,7 @@ int back_prevalence_limit(double *p, dou fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp); } /* Age */ /* was end of cptcod */ + /*fprintf(ficresplb,"\n");*/ /* Seems to be necessary for gnuplot only if two result lines and no covariate. */ } /* end of any combination */ } /* end of nres */ /* hBijx(p, bage, fage); */ @@ -9367,7 +9614,7 @@ int hPijx(double *p, int bage, int fage) /* k=k+1; */ for(nres=1; nres <= nresult; nres++) /* For each resultline */ for(k=1; k<=i1;k++){ - if(TKresult[nres]!= k) + if(i1 != 1 && TKresult[nres]!= k) continue; fprintf(ficrespij,"\n#****** "); for(j=1;j<=cptcoveff;j++) @@ -9439,14 +9686,14 @@ 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 "); + fprintf(ficrespijb,"#****** h Bij x Back probability to be in state i at age x-h being in j at x: B1j+B2j+...=1 "); i1= pow(2,cptcoveff); /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */ /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */ /* k=k+1; */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ - if(TKresult[nres]!= k) + if(i1 != 1 && TKresult[nres]!= k) continue; fprintf(ficrespijb,"\n#****** "); for(j=1;j<=cptcoveff;j++) @@ -9473,7 +9720,7 @@ int hPijx(double *p, int bage, int fage) /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */ hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k); /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */ - fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j="); + fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j="); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespijb," %1d-%1d",i,j); @@ -9563,7 +9810,8 @@ int main(int argc, char *argv[]) double **prlim; double **bprlim; double ***param; /* Matrix of parameters */ - double *p; + double ***paramstart; /* Matrix of starting parameter values */ + double *p, *pstart; /* p=param[1][1] pstart is for starting values guessed by freqsummary */ double **matcov; /* Matrix of covariance */ double **hess; /* Hessian matrix */ double ***delti3; /* Scale */ @@ -9874,6 +10122,12 @@ int main(int argc, char *argv[]) delti=delti3[1][1]; /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ if(mle==-1){ /* Print a wizard for help writing covariance matrix */ +/* We could also provide initial parameters values giving by simple logistic regression + * only one way, that is without matrix product. We will have nlstate maximizations */ + /* 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, \ + freqsummary(fileres, p, pstart, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ firstpass, lastpass, stepm, weightopt, model); fprintf(fichtm,"\n"); @@ -10408,9 +10664,9 @@ Interval (in months) between two waves: /* For mortality only */ if (mle==-3){ ximort=matrix(1,NDIM,1,NDIM); - for(i=1;i<=NDIM;i++) - for(j=1;j<=NDIM;j++) - ximort[i][j]=0.; + for(i=1;i<=NDIM;i++) + for(j=1;j<=NDIM;j++) + ximort[i][j]=0.; /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ cens=ivector(1,n); ageexmed=vector(1,n); @@ -10646,6 +10902,10 @@ Please run with mle=-1 to get a correct printf("\n"); if(mle>=1){ /* Could be 1 or 2, Real Maximization */ /* mlikeli uses func not funcone */ + /* for(i=1;i