--- imach/src/imach.c 2014/09/25 11:43:39 1.162 +++ imach/src/imach.c 2024/04/30 10:59:22 1.360 @@ -1,689 +1,1554 @@ -/* $Id: imach.c,v 1.162 2014/09/25 11:43:39 brouard Exp $ +/* $Id: imach.c,v 1.360 2024/04/30 10:59:22 brouard Exp $ $State: Exp $ $Log: imach.c,v $ - Revision 1.162 2014/09/25 11:43:39 brouard - Summary: temporary backup 0.99! + Revision 1.360 2024/04/30 10:59:22 brouard + Summary: Version 0.99s4 and estimation of std of e.j/e.. - Revision 1.1 2014/09/16 11:06:58 brouard - Summary: With some code (wrong) for nlopt + Revision 1.359 2024/04/24 21:21:17 brouard + Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes - Author: + Revision 1.6 2024/04/24 21:10:29 brouard + Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes - Revision 1.161 2014/09/15 20:41:41 brouard - Summary: Problem with macro SQR on Intel compiler + Revision 1.5 2023/10/09 09:10:01 brouard + Summary: trying to reconsider - Revision 1.160 2014/09/02 09:24:05 brouard + Revision 1.4 2023/06/22 12:50:51 brouard + Summary: stil on going + + Revision 1.3 2023/06/22 11:28:07 brouard *** empty log message *** - Revision 1.159 2014/09/01 10:34:10 brouard - Summary: WIN32 - Author: Brouard + Revision 1.2 2023/06/22 11:22:40 brouard + Summary: with svd but not working yet - Revision 1.158 2014/08/27 17:11:51 brouard + Revision 1.353 2023/05/08 18:48:22 brouard *** empty log message *** - Revision 1.157 2014/08/27 16:26:55 brouard - Summary: Preparing windows Visual studio version - Author: Brouard + Revision 1.352 2023/04/29 10:46:21 brouard + *** empty log message *** - In order to compile on Visual studio, time.h is now correct and time_t - and tm struct should be used. difftime should be used but sometimes I - just make the differences in raw time format (time(&now). - Trying to suppress #ifdef LINUX - Add xdg-open for __linux in order to open default browser. + Revision 1.351 2023/04/29 10:43:47 brouard + Summary: 099r45 - Revision 1.156 2014/08/25 20:10:10 brouard + Revision 1.350 2023/04/24 11:38:06 brouard *** empty log message *** - Revision 1.155 2014/08/25 18:32:34 brouard - Summary: New compile, minor changes - Author: Brouard + Revision 1.349 2023/01/31 09:19:37 brouard + Summary: Improvements in models with age*Vn*Vm - Revision 1.154 2014/06/20 17:32:08 brouard - Summary: Outputs now all graphs of convergence to period prevalence + Revision 1.347 2022/09/18 14:36:44 brouard + Summary: version 0.99r42 - Revision 1.153 2014/06/20 16:45:46 brouard - Summary: If 3 live state, convergence to period prevalence on same graph - Author: Brouard + Revision 1.346 2022/09/16 13:52:36 brouard + * src/imach.c (Module): 0.99r41 Was an error when product of timevarying and fixed. Using FixedV[of name] now. Thank you Feinuo - Revision 1.152 2014/06/18 17:54:09 brouard - Summary: open browser, use gnuplot on same dir than imach if not found in the path + Revision 1.345 2022/09/16 13:40:11 brouard + Summary: Version 0.99r41 - Revision 1.151 2014/06/18 16:43:30 brouard - *** empty log message *** + * imach.c (Module): 0.99r41 Was an error when product of timevarying and fixed. Using FixedV[of name] now. Thank you Feinuo - Revision 1.150 2014/06/18 16:42:35 brouard - Summary: If gnuplot is not in the path try on same directory than imach binary (OSX) - Author: brouard + Revision 1.344 2022/09/14 19:33:30 brouard + Summary: version 0.99r40 - Revision 1.149 2014/06/18 15:51:14 brouard - Summary: Some fixes in parameter files errors - Author: Nicolas Brouard + * imach.c (Module): Fixing names of variables in T_ (thanks to Feinuo) - Revision 1.148 2014/06/17 17:38:48 brouard - Summary: Nothing new - Author: Brouard + Revision 1.343 2022/09/14 14:22:16 brouard + Summary: version 0.99r39 - Just a new packaging for OS/X version 0.98nS + * imach.c (Module): Version 0.99r39 with colored dummy covariates + (fixed or time varying), using new last columns of + ILK_parameter.txt file. - Revision 1.147 2014/06/16 10:33:11 brouard - *** empty log message *** + Revision 1.342 2022/09/11 19:54:09 brouard + Summary: 0.99r38 - Revision 1.146 2014/06/16 10:20:28 brouard - Summary: Merge - Author: Brouard + * imach.c (Module): Adding timevarying products of any kinds, + should work before shifting cotvar from ncovcol+nqv columns in + order to have a correspondance between the column of cotvar and + the id of column. + (Module): Some cleaning and adding covariates in ILK.txt - Merge, before building revised version. + Revision 1.341 2022/09/11 07:58:42 brouard + Summary: Version 0.99r38 - Revision 1.145 2014/06/10 21:23:15 brouard - Summary: Debugging with valgrind - Author: Nicolas Brouard + After adding change in cotvar. - Lot of changes in order to output the results with some covariates - After the Edimburgh REVES conference 2014, it seems mandatory to - improve the code. - No more memory valgrind error but a lot has to be done in order to - continue the work of splitting the code into subroutines. - Also, decodemodel has been improved. Tricode is still not - optimal. nbcode should be improved. Documentation has been added in - the source code. + Revision 1.340 2022/09/11 07:53:11 brouard + Summary: Version imach 0.99r37 - Revision 1.143 2014/01/26 09:45:38 brouard - Summary: Version 0.98nR (to be improved, but gives same optimization results as 0.98k. Nice, promising + * imach.c (Module): Adding timevarying products of any kinds, + should work before shifting cotvar from ncovcol+nqv columns in + order to have a correspondance between the column of cotvar and + the id of column. - * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... - (Module): Version 0.98nR Running ok, but output format still only works for three covariates. + Revision 1.339 2022/09/09 17:55:22 brouard + Summary: version 0.99r37 - Revision 1.142 2014/01/26 03:57:36 brouard - Summary: gnuplot changed plot w l 1 has to be changed to plot w l lt 2 + * imach.c (Module): Many improvements for fixing products of fixed + timevarying as well as fixed * fixed, and test with quantitative + covariate. - * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... + Revision 1.338 2022/09/04 17:40:33 brouard + Summary: 0.99r36 - Revision 1.141 2014/01/26 02:42:01 brouard - * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... + * imach.c (Module): Now the easy runs i.e. without result or + model=1+age only did not work. The defautl combination should be 1 + and not 0 because everything hasn't been tranformed yet. - Revision 1.140 2011/09/02 10:37:54 brouard - Summary: times.h is ok with mingw32 now. + Revision 1.337 2022/09/02 14:26:02 brouard + Summary: version 0.99r35 - Revision 1.139 2010/06/14 07:50:17 brouard - After the theft of my laptop, I probably lost some lines of codes which were not uploaded to the CVS tree. - I remember having already fixed agemin agemax which are pointers now but not cvs saved. + * src/imach.c: Version 0.99r35 because it outputs same results with + 1+age+V1+V1*age for females and 1+age for females only + (education=1 noweight) - Revision 1.138 2010/04/30 18:19:40 brouard + Revision 1.336 2022/08/31 09:52:36 brouard *** empty log message *** - Revision 1.137 2010/04/29 18:11:38 brouard - (Module): Checking covariates for more complex models - than V1+V2. A lot of change to be done. Unstable. + Revision 1.335 2022/08/31 08:23:16 brouard + Summary: improvements... - Revision 1.136 2010/04/26 20:30:53 brouard - (Module): merging some libgsl code. Fixing computation - of likelione (using inter/intrapolation if mle = 0) in order to - get same likelihood as if mle=1. - Some cleaning of code and comments added. + Revision 1.334 2022/08/25 09:08:41 brouard + Summary: In progress for quantitative - Revision 1.135 2009/10/29 15:33:14 brouard - (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code. + Revision 1.333 2022/08/21 09:10:30 brouard + * src/imach.c (Module): Version 0.99r33 A lot of changes in + reassigning covariates: my first idea was that people will always + use the first covariate V1 into the model but in fact they are + producing data with many covariates and can use an equation model + with some of the covariate; it means that in a model V2+V3 instead + of codtabm(k,Tvaraff[j]) which calculates for combination k, for + three covariates (V1, V2, V3) the value of Tvaraff[j], but in fact + the equation model is restricted to two variables only (V2, V3) + and the combination for V2 should be codtabm(k,1) instead of + (codtabm(k,2), and the code should be + codtabm(k,TnsdVar[Tvaraff[j]]. Many many changes have been + made. All of these should be simplified once a day like we did in + hpxij() for example by using precov[nres] which is computed in + decoderesult for each nres of each resultline. Loop should be done + on the equation model globally by distinguishing only product with + age (which are changing with age) and no more on type of + covariates, single dummies, single covariates. + + Revision 1.332 2022/08/21 09:06:25 brouard + Summary: Version 0.99r33 + + * src/imach.c (Module): Version 0.99r33 A lot of changes in + reassigning covariates: my first idea was that people will always + use the first covariate V1 into the model but in fact they are + producing data with many covariates and can use an equation model + with some of the covariate; it means that in a model V2+V3 instead + of codtabm(k,Tvaraff[j]) which calculates for combination k, for + three covariates (V1, V2, V3) the value of Tvaraff[j], but in fact + the equation model is restricted to two variables only (V2, V3) + and the combination for V2 should be codtabm(k,1) instead of + (codtabm(k,2), and the code should be + codtabm(k,TnsdVar[Tvaraff[j]]. Many many changes have been + made. All of these should be simplified once a day like we did in + hpxij() for example by using precov[nres] which is computed in + decoderesult for each nres of each resultline. Loop should be done + on the equation model globally by distinguishing only product with + age (which are changing with age) and no more on type of + covariates, single dummies, single covariates. - Revision 1.134 2009/10/29 13:18:53 brouard - (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code. + Revision 1.331 2022/08/07 05:40:09 brouard + *** empty log message *** - Revision 1.133 2009/07/06 10:21:25 brouard - just nforces + Revision 1.330 2022/08/06 07:18:25 brouard + Summary: last 0.99r31 - Revision 1.132 2009/07/06 08:22:05 brouard - Many tings + * imach.c (Module): Version of imach using partly decoderesult to rebuild xpxij function - Revision 1.131 2009/06/20 16:22:47 brouard - Some dimensions resccaled + Revision 1.329 2022/08/03 17:29:54 brouard + * imach.c (Module): Many errors in graphs fixed with Vn*age covariates. - Revision 1.130 2009/05/26 06:44:34 brouard - (Module): Max Covariate is now set to 20 instead of 8. A - lot of cleaning with variables initialized to 0. Trying to make - V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better. + Revision 1.328 2022/07/27 17:40:48 brouard + Summary: valgrind bug fixed by initializing to zero DummyV as well as Tage - Revision 1.129 2007/08/31 13:49:27 lievre - Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting + Revision 1.327 2022/07/27 14:47:35 brouard + Summary: Still a problem for one-step probabilities in case of quantitative variables - Revision 1.128 2006/06/30 13:02:05 brouard - (Module): Clarifications on computing e.j + Revision 1.326 2022/07/26 17:33:55 brouard + Summary: some test with nres=1 - Revision 1.127 2006/04/28 18:11:50 brouard - (Module): Yes the sum of survivors was wrong since - imach-114 because nhstepm was no more computed in the age - loop. Now we define nhstepma in the age loop. - (Module): In order to speed up (in case of numerous covariates) we - compute health expectancies (without variances) in a first step - and then all the health expectancies with variances or standard - deviation (needs data from the Hessian matrices) which slows the - computation. - In the future we should be able to stop the program is only health - expectancies and graph are needed without standard deviations. + Revision 1.325 2022/07/25 14:27:23 brouard + Summary: r30 - Revision 1.126 2006/04/28 17:23:28 brouard - (Module): Yes the sum of survivors was wrong since - imach-114 because nhstepm was no more computed in the age - loop. Now we define nhstepma in the age loop. - Version 0.98h + * imach.c (Module): Error cptcovn instead of nsd in bmij (was + coredumped, revealed by Feiuno, thank you. - Revision 1.125 2006/04/04 15:20:31 lievre - Errors in calculation of health expectancies. Age was not initialized. - Forecasting file added. + Revision 1.324 2022/07/23 17:44:26 brouard + *** empty log message *** - Revision 1.124 2006/03/22 17:13:53 lievre - Parameters are printed with %lf instead of %f (more numbers after the comma). - The log-likelihood is printed in the log file + Revision 1.323 2022/07/22 12:30:08 brouard + * imach.c (Module): Output of Wald test in the htm file and not only in the log. - Revision 1.123 2006/03/20 10:52:43 brouard - * imach.c (Module):
\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); + } + /* for (z1=1; z1<= nqtveff; z1++) { */ + /* for(m=1;m<=lastpass;m++){ */ + /* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f
\n", z1, m, meanqt[m][z1]); */ + /* } */ + /* } */ + + fprintf(ficresphtm,"
Age | Prev(%d) | N(%d) | N | ",i,i); } - for(i=1; i<=nlstate;i++) - fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); - fprintf(ficresp, "\n"); + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); + fprintf(ficresphtm, "\n"); - for(i=iagemin; i <= iagemax+3; i++){ - if(i==iagemax+3){ + /* Header of frequency table by age */ + fprintf(ficresphtmfr,"
---|
Age | "); + for(s2=-1; s2 <=nlstate+ndeath; s2++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(s2!=0 && m!=0) + fprintf(ficresphtmfr,"%d%d | ",s2,m); + } + } + fprintf(ficresphtmfr, "\n"); + + /* For each age */ + for(iage=iagemin; iage <= iagemax+3; iage++){ + fprintf(ficresphtm,"||||||||
---|---|---|---|---|---|---|---|---|---|
0 | "); + }else if(iage==iagemax+2){ + fprintf(ficlog,"0"); + fprintf(ficresphtmfr,"|||||||||
Unknown | "); + }else if(iage==iagemax+3){ fprintf(ficlog,"Total"); + fprintf(ficresphtmfr,"|||||||||
Total | "); }else{ if(first==1){ first=0; printf("See log file for details...\n"); } - fprintf(ficlog,"Age %d", i); + fprintf(ficresphtmfr,"|||||||||
%d | ",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][i]; + for(s1=1; s1 <=nlstate ; s1++){ + for(m=-1, pp[s1]=0; m <=nlstate+ndeath ; m++) + pp[s1] += freq[s1][m][iage]; } - for(jk=1; jk <=nlstate ; jk++){ + for(s1=1; s1 <=nlstate ; s1++){ for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][i]; - if(pp[jk]>=1.e-10){ + pos += freq[s1][m][iage]; + if(pp[s1]>=1.e-10){ if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + printf(" %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); } - fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); }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); + printf(" %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); + fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); } } - - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1,pos=0,posprop=0; jk <=nlstate ; jk++){ - pos += pp[jk]; - posprop += prop[jk][i]; + + for(s1=1; s1 <=nlstate ; s1++){ + /* posprop[s1]=0; */ + for(m=0, pp[s1]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ + pp[s1] += freq[s1][m][iage]; + } /* pp[s1] is the total number of transitions starting from state s1 and any ending status until this age */ + + for(s1=1,pos=0, pospropta=0.; s1 <=nlstate ; s1++){ + pos += pp[s1]; /* pos is the total number of transitions until this age */ + posprop[s1] += prop[s1][iage]; /* prop is the number of transitions from a live state + from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pospropta += prop[s1][iage]; /* prop is the number of transitions from a live state + from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + } + + /* Writing ficresp */ + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + if( iage <= iagemax){ + fprintf(ficresp," %d",iage); + } + }else if( nj==2){ + if( iage <= iagemax){ + fprintf(ficresp," %d",iage); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + } } - for(jk=1; jk <=nlstate ; jk++){ + for(s1=1; s1 <=nlstate ; s1++){ 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); + printf(" %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos); + fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/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); + printf(" %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); + fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); } - if( i <= iagemax){ + if( iage <= iagemax){ if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); - /*probs[i][jk][j1]= pp[jk]/pos;*/ - /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + }else if( nj==2){ + fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + } + fprintf(ficresphtm,"%d | %.5f | %.0f | %.0f | ",iage,prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + /*probs[iage][s1][j1]= pp[s1]/pos;*/ + /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/ + } else{ + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); + fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f | ",iage, prop[s1][iage],pospropta); } - else - fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); } - } - - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) - if(freq[jk][m][i] !=0 ) { - if(first==1) - printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]); + pospropt[s1] +=posprop[s1]; + } /* end loop s1 */ + /* pospropt=0.; */ + for(s1=-1; s1 <=nlstate+ndeath; s1++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[s1][m][iage] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); + } + /* printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); */ + fprintf(ficlog," %d%d=%.0f",s1,m,freq[s1][m][iage]); } - if(i <= iagemax) - fprintf(ficresp,"\n"); + if(s1!=0 && m!=0) + fprintf(ficresphtmfr,"%.0f | ",freq[s1][m][iage]); + } + } /* end loop s1 */ + posproptt=0.; + for(s1=1; s1 <=nlstate; s1++){ + posproptt += pospropt[s1]; + } + fprintf(ficresphtmfr,"
Tot | "); + for(s1=1; s1 <=nlstate ; s1++){ + if(posproptt < 1.e-5){ + fprintf(ficresphtm,"Nanq | %.0f | %.0f | ",pospropt[s1],posproptt); + }else{ + fprintf(ficresphtm,"%.5f | %.0f | %.0f | ",pospropt[s1]/posproptt,pospropt[s1],posproptt); + } } - /*}*/ - } + fprintf(ficresphtm,"
This combination (%d) is not valid and no result will be produced
",j1); + fprintf(ficresphtmfr,"\nThis combination (%d) is not valid and no result will be produced
",j1); + fprintf(ficlog,"# This combination (%d) is not valid and no result will be produced\n",j1); + printf("# This combination (%d) is not valid and no result will be produced\n",j1); + invalidvarcomb[j1]=1; + }else{ + fprintf(ficresphtm,"\nThis combination (%d) is valid and result will be produced (or no resultline).
",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,s1=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + for(jj=1; jj <=ncovmodel; jj++){ /* For counting s1 */ + 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[s1],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[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + pstart[s1]= 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, s1=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,s1,j1,jj, iage, y[iage]); */ + } + /* Some are not finite, but linreg will ignore these ages */ + no=0; + linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */ + pstart[s1]=b; + pstart[s1-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[s1]= 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("s1=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",s1,i,k,s1,p[s1],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]); */ + s1++; + } /* end jj */ + } /* end k!= i */ + } /* end k */ + } /* end i, s1 */ + } /* 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,s1=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[s1]=p[s1]; /* Setting pstart to p values by default */ + if(jj==1){ /* Age has to be done */ + pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); + printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],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[s1],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]); */ + s1++; + } + printf("\n"); + fprintf(ficlog,"\n"); + } + } + } /* end of state i */ + printf("#Freqsummary\n"); + fprintf(ficlog,"\n"); + for(s1=-1; s1 <=nlstate+ndeath; s1++){ + for(s2=-1; s2 <=nlstate+ndeath; s2++){ + /* param[i]|j][k]= freq[s1][s2][iagemax+3] */ + printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); + fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); + /* if(freq[s1][s2][iage] !=0 ) { /\* minimizing output *\/ */ + /* printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ + /* fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ + /* } */ + } + } /* end loop s1 */ + + 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"); + } + } + } + } /* end mle=-2 */ dateintmean=dateintsum/k2cpt; - + date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); + fclose(ficresp); - free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); + fclose(ficresphtm); + fclose(ficresphtmfr); + free_vector(idq,1,nqfveff); + free_vector(meanq,1,nqfveff); + free_vector(stdq,1,nqfveff); + free_matrix(meanqt,1,lastpass,1,nqtveff); + 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+4+AGEMARGE); free_vector(pp,1,nlstate); - free_matrix(prop,1,nlstate,iagemin, iagemax+3); - /* End of Freq */ + /* End of freqsummary */ +} + +/* Simple linear regression */ +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 = 0.0; /* 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/(double)(ne-2)/(sumx2 - sumx * sumx /(double)ne)); + *sa= *sb * sqrt(sumx2/ne); + + return 0; } /************ Prevalence ********************/ @@ -2650,10 +7708,13 @@ void prevalence(double ***probs, double We still use firstpass and lastpass as another selection. */ - int i, m, jk, k1, i1, j1, bool, z1,j; - double ***freq; /* Frequencies */ - double *pp, **prop; - double pos,posprop; + int i, m, jk, j1, bool, z1,j, iv; + int mi; /* Effective wave */ + int iage; + double agebegin; /*, ageend;*/ + + double **prop; + double posprop; double y2; /* in fractional years */ int iagemin, iagemax; int first; /** to stop verbosity which is redirected to log file */ @@ -2661,111 +7722,206 @@ void prevalence(double ***probs, double iagemin= (int) agemin; iagemax= (int) agemax; /*pp=vector(1,nlstate);*/ - prop=matrix(1,nlstate,iagemin,iagemax+3); + prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ j1=0; /*j=cptcoveff;*/ if (cptcovn<1) {j=1;ncodemax[1]=1;} - first=1; - for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ - /*for(i1=1; i1<=ncodemax[k1];i1++){ - j1++;*/ - - for (i=1; i<=nlstate; i++) - for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0.0; - - for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + first=0; + for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of simple dummy covariates */ + for (i=1; i<=nlstate; i++) + 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)]); */ + fprintf(ficlog,"Prevalence combination of varying and fixed dummies %d\n",j1); + + for (i=1; i<=imx; i++) { /* Each individual */ + bool=1; + /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ + for(mi=1; mi
Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\
+ : \
+%s_%d%1d%1d-%1d%1d.svg, ",k1,l1,k2,l2,\
+ subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \
+ subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+ fprintf(fichtmcov,"\n
",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+ fprintf(fichtmcov,"\n
Correlation at age %d (%.3f),",(int) age, c12);
+ fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+ fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
+ fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
+ fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \
+ mu1,std,v11,sqrt(fabs(lc1)),v12,sqrt(fabs(lc2)), \
+ mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */
+ }else{
+ first=0;
+ fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
+ fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
+ fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
+ fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \
+ mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \
+ mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2)));
+ }/* if first */
+ } /* age mod 5 */
+ } /* end loop age */
+ fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+ first=1;
+ } /*l12 */
+ } /* k12 */
+ } /*l1 */
+ }/* k1 */
+ } /* loop on combination of covariates j1 */
+ } /* loop on nres */
+ free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
+ free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
+ free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
+ free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar);
+ free_vector(xp,1,npar);
+ fclose(ficresprob);
+ fclose(ficresprobcov);
+ fclose(ficresprobcor);
+ fflush(ficgp);
+ fflush(fichtmcov);
+ }
/******************* Printing html file ***********/
-void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
+void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
int lastpass, int stepm, int weightopt, char model[],\
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
- int popforecast, int estepm ,\
- double jprev1, double mprev1,double anprev1, \
- double jprev2, double mprev2,double anprev2){
- int jj1, k1, i1, cpt;
-
+ int popforecast, int mobilav, int prevfcast, int mobilavproj, int prevbcast, int estepm , \
+ double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \
+ double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){
+ int jj1, k1, cpt, nres;
+ /* In fact some results are already printed in fichtm which is open */
fprintf(fichtm,"
"); + + jj1=0; + + fprintf(fichtm," \n
"); + jj1=0; - m=pow(2,cptcoveff); - if (cptcovn < 1) {m=1;ncodemax[1]=1;} + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ + /* k1=nres; */ + k1=TKresult[nres]; + if(TKresult[nres]==0) k1=1; /* To be checked for noresult */ + /* for(k1=1; k1<=m;k1++){ /\* For each combination of covariate *\/ */ + /* if(m != 1 && TKresult[nres]!= k1) */ + /* continue; */ - jj1=0; - for(k1=1; k1<=m;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ + /* for(i1=1; i1<=ncodemax[k1];i1++){ */ jj1++; if (cptcovn > 0) { + fprintf(fichtm,"\n
");
+
fprintf(fichtm,"
************ Results for covariates");
- for (cpt=1; cpt<=cptcoveff;cpt++)
- fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
- fprintf(fichtm," ************\n
");
+ for (cpt=1; cpt<=cptcovs;cpt++){
+ fprintf(fichtm," V%d=%lg ",Tvresult[nres][cpt],TinvDoQresult[nres][Tvresult[nres][cpt]]);
+ printf(" V%d=%lg ",Tvresult[nres][cpt],TinvDoQresult[nres][Tvresult[nres][cpt]]);
+ /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+ /* printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); */
+ }
+ /* if(nqfveff+nqtveff 0) */ /* Test to be done */
+ fprintf(fichtm," (model=1+age+%s) ************\n
",model);
+ if(invalidvarcomb[k1]){
+ fprintf(fichtm,"\nCombination (%d) ignored because no cases
\n",k1);
+ printf("\nCombination (%d) ignored because no cases \n",k1);
+ continue;
+ }
}
+ /* aij, bij */
+ fprintf(fichtm,"
- Logit model (yours is: logit(pij)=log(pij/pii)= aij+ bij age+%s) as a function of age: %s_%d-1-%d.svg
\
+",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Pij */
- fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d_1.png
\
-",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ fprintf(fichtm,"
\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s_%d-2-%d.svg
\
+",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Quasi-incidences */
- fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
- before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d_2.png
\
-",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
- /* Period (stable) prevalence in each health state */
+ fprintf(fichtm,"
\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \
+ incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \
+divided by h: hPij/h : %s_%d-3-%d.svg
\
+",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
+ /* Survival functions (period) in state j */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. %s_%d-%d-%d.svg
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
+ fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+ }
+ /* State specific survival functions (period) */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Survival functions in state %d and in any other live state (total).\
+ And probability to be observed in various states (up to %d) being in state %d at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. \
+ %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
+ fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+ }
+ /* Period (forward 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 for a person being in state (1 to %d) at different ages, to be alive in state %d some years after. %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
+ fprintf(fichtm,"" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+ }
+ if(prevbcast==1){
+ /* Backward prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was 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);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJB_"),subdirf2(optionfilefiname,"PIJB_"));
+ fprintf(fichtm,"" ,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
+ }
+ }
+ if(prevfcast==1){
+ /* Projection of prevalence up to period (forward stable) prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg", dateprev1, dateprev2, mobilavproj, dateprojd, dateprojf, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"F_"),subdirf2(optionfilefiname,"F_"));
+ fprintf(fichtm,"",
+ subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+ }
+ }
+ if(prevbcast==1){
+ /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"
- Convergence from each state (1 to %d) to period (stable) prevalence in state %d %s%d_%d.png
\
-",nlstate, cpt, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+ fprintf(fichtm,"
\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
+ from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d). Randomness in cross-sectional prevalence is not taken into \
+ account but can visually be appreciated. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
+with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));
+ fprintf(fichtm," ", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
}
+ }
+
for(cpt=1; cpt<=nlstate;cpt++) {
- fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : %s%d%d.png
\
-",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+ fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d-%d-%d.svg",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"E_"),subdirf2(optionfilefiname,"E_"));
+ fprintf(fichtm,"", subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres );
}
- } /* end i1 */
- }/* End k1 */
- fprintf(fichtm,"
");
-
- m=pow(2,cptcoveff);
- if (cptcovn < 1) {m=1;ncodemax[1]=1;}
-
- jj1=0;
- for(k1=1; k1<=m;k1++){
- for(i1=1; i1<=ncodemax[k1];i1++){
+/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=1+age+%s (instead of .)
\n",popforecast, stepm, model); */
+ fflush(fichtm);
+
+ m=pow(2,cptcoveff);
+ if (cptcovn < 1) {m=1;ncodemax[1]=1;}
+
+ fprintf(fichtm,"
"); + + jj1=0; + + fprintf(fichtm," \n
");
+
fprintf(fichtm,"
************ Results for covariates");
- for (cpt=1; cpt<=cptcoveff;cpt++)
- fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
- fprintf(fichtm," ************\n
");
- }
+ for (cpt=1; cpt<=cptcovs;cpt++){ /**< cptcoveff number of variables */
+ fprintf(fichtm," V%d=%lg ",Tvresult[nres][cpt],TinvDoQresult[nres][Tvresult[nres][cpt]]);
+ printf(" V%d=%lg ",Tvresult[nres][cpt],TinvDoQresult[nres][Tvresult[nres][cpt]]);
+ /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+ }
+
+ fprintf(fichtm," (model=1+age+%s) ************\n
",model);
+
+ if(invalidvarcomb[k1]){
+ fprintf(fichtm,"\nCombination (%d) ignored because no cases
\n",k1);
+ continue;
+ }
+ } /* If cptcovn >0 */
for(cpt=1; cpt<=nlstate;cpt++) {
- fprintf(fichtm,"
- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.png
\
-",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
+ fprintf(fichtm,"\n
- Observed (cross-sectional with mov_average=%d) and period (incidence based) \
+prevalence (with 95%% confidence interval) in state (%d): %s_%d-%d-%d.svg",mobilav,cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s)\n
",subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
+ fprintf(fichtm,"",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
}
fprintf(fichtm,"\n
- Total life expectancy by age and \
-health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
+health expectancies in each live state (1 to %d) with confidence intervals \
+on left y-scale as well as proportions of time spent in each live state \
+(with confidence intervals) on right y-scale 0 to 100%%.\
+ If popbased=1 the smooth (due to the model) \
true period expectancies (those weighted with period prevalences are also\
drawn in addition to the population based expectancies computed using\
- observed and cahotic prevalences: %s%d.png
\
-",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
- } /* end i1 */
- }/* End k1 */
- fprintf(fichtm,"
Model= | 1 | + age | "); + if(nagesqr==1){ + printf(" + age*age "); + fprintf(ficres," + age*age "); + fprintf(ficlog," + age*age "); + fprintf(fichtm, "+ age*age | "); + } + for(j=1;j <=ncovmodel-2;j++){ + if(Typevar[j]==0) { + printf(" + V%d ",Tvar[j]); + fprintf(ficres," + V%d ",Tvar[j]); + fprintf(ficlog," + V%d ",Tvar[j]); + fprintf(fichtm, "+ V%d | ",Tvar[j]); + }else if(Typevar[j]==1) { + printf(" + V%d*age ",Tvar[j]); + fprintf(ficres," + V%d*age ",Tvar[j]); + fprintf(ficlog," + V%d*age ",Tvar[j]); + fprintf(fichtm, "+ V%d*age | ",Tvar[j]); + }else if(Typevar[j]==2) { + printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficres," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(fichtm, "+ V%d*V%d | ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + }else if(Typevar[j]==3) { /* TO VERIFY */ + printf(" + V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficres," + V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficlog," + V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(fichtm, "+ V%d*V%d*age | ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + } + } + printf("\n"); + fprintf(ficres,"\n"); + fprintf(ficlog,"\n"); + fprintf(fichtm, "
---|---|---|---|---|---|---|---|
%1d%1d | ",i,k); for(j=1; j <=ncovmodel; j++){ - printf("%lf ",p[jk]); - fprintf(ficlog,"%lf ",p[jk]); - fprintf(ficres,"%lf ",p[jk]); + printf("%12.7f ",p[jk]); + fprintf(ficlog,"%12.7f ",p[jk]); + fprintf(ficres,"%12.7f ",p[jk]); + fprintf(fichtm, "%12.7f | ",p[jk]); jk++; } printf("\n"); fprintf(ficlog,"\n"); fprintf(ficres,"\n"); + fprintf(fichtm, "
The Wald test results are output only if the maximimzation of the Likelihood is performed (mle=1)\nParameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities and its graphs'.\n",optionfilehtmcov); + fprintf(fichtm,"\n
Model= | 1 | + age | "); + if(nagesqr==1){ + printf(" + age*age "); + fprintf(ficres," + age*age "); + fprintf(ficlog," + age*age "); + fprintf(fichtm, "+ age*age | "); + } + for(j=1;j <=ncovmodel-2;j++){ + if(Typevar[j]==0) { + printf(" + V%d ",Tvar[j]); + fprintf(fichtm, "+ V%d | ",Tvar[j]); + }else if(Typevar[j]==1) { + printf(" + V%d*age ",Tvar[j]); + fprintf(fichtm, "+ V%d*age | ",Tvar[j]); + }else if(Typevar[j]==2) { + fprintf(fichtm, "+ V%d*V%d | ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + }else if(Typevar[j]==3) { /* TO VERIFY */ + fprintf(fichtm, "+ V%d*V%d*age | ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + } + } + fprintf(fichtm, "
---|---|---|---|---|---|---|---|
%1d%1d | ",i,k); + for(j=1; j <=ncovmodel; j++){ + wald=p[jk]/sqrt(matcov[jk][jk]); + printf("%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); + fprintf(ficlog,"%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); + if(fabs(wald) > 1.96){ + fprintf(fichtm, "%12.7f (%12.7f)",p[jk],sqrt(matcov[jk][jk])); + }else{ + fprintf(fichtm, " | %12.7f (%12.7f)",p[jk],sqrt(matcov[jk][jk])); + } + fprintf(fichtm,"W=%8.3f",wald); + fprintf(fichtm,"[%12.7f;%12.7f] | ", p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); + jk++; + } + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(fichtm, "