/* $Id$
$State$
$Log$
+ Revision 1.225 2016/07/12 08:40:03 brouard
+ Summary: saving but not running
+
Revision 1.224 2016/07/01 13:16:01 brouard
Summary: Fixes
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
nhstepm*hstepm matrices. Returns p3mat[i][j][h] after calling
p3mat[i][j][h]=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\
- 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+ 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+
+Important routines
+
+- func (or funcone), computes logit (pij) distinguishing
+ o fixed variables (single or product dummies or quantitative);
+ o varying variables by:
+ (1) wave (single, product dummies, quantitative),
+ (2) by age (can be month) age (done), age*age (done), age*Vn where Vn can be:
+ % fixed dummy (treated) or quantitative (not done because time-consuming);
+ % varying dummy (not done) or quantitative (not done);
+- Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)
+ and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually.
+- printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables
+ o There are 2*cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if
+ race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.
+
+
Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
Institut national d'études démographiques, Paris.
This software have been partly granted by Euro-REVES, a concerted action
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <ctype.h>
#ifdef _WIN32
#include <io.h>
double ***cotqvar; /* Time varying quantitative covariate itqv */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
-int *Typevar; /**< 1 for qualitative fixed, 2 for quantitative fixed, 3 for qualitive varying, 4 for quanti varying*/
+int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */
+int *Fixed; /** Fixed[Tvar[k]] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
+int *Dummy; /** Dummy[Tvar[k]] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
int *Tage;
int *Ndum; /** Freq of modality (tricode */
/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
/*************** log-likelihood *************/
double func( double *x)
{
- int i, ii, j, k, mi, d, kk;
- int ioffset=0;
- double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
- double **out;
- double sw; /* Sum of weights */
- double lli; /* Individual log likelihood */
- int s1, s2;
- int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quatitative time varying covariate */
- double bbh, survp;
- long ipmx;
- double agexact;
- /*extern weight */
- /* We are differentiating ll according to initial status */
- /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
- /*for(i=1;i<imx;i++)
- printf(" %d\n",s[4][i]);
- */
+ int i, ii, j, k, mi, d, kk;
+ int ioffset=0;
+ double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
+ double **out;
+ double sw; /* Sum of weights */
+ double lli; /* Individual log likelihood */
+ int s1, s2;
+ int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quatitative time varying covariate */
+ double bbh, survp;
+ long ipmx;
+ double agexact;
+ /*extern weight */
+ /* We are differentiating ll according to initial status */
+ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
+ /*for(i=1;i<imx;i++)
+ printf(" %d\n",s[4][i]);
+ */
- ++countcallfunc;
+ ++countcallfunc;
- cov[1]=1.;
+ cov[1]=1.;
- for(k=1; k<=nlstate; k++) ll[k]=0.;
+ for(k=1; k<=nlstate; k++) ll[k]=0.;
ioffset=0;
- if(mle==1){
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- /* Computes the values of the ncovmodel covariates of the model
- depending if the covariates are fixed or varying (age dependent) and stores them in cov[]
- Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
- to be observed in j being in i according to the model.
- */
- ioffset=2+nagesqr+cptcovage;
- /* for (k=1; k<=cptcovn;k++){ /\* Simple and product covariates without age* products *\/ */
- for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */
- cov[++ioffset]=covar[Tvar[k]][i];
- }
- for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitatives and Fixed covariates */
- cov[++ioffset]=coqvar[iqv][i];
- }
+ if(mle==1){
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ /* Computes the values of the ncovmodel covariates of the model
+ depending if the covariates are fixed or varying (age dependent) and stores them in cov[]
+ Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
+ to be observed in j being in i according to the model.
+ */
+ ioffset=2+nagesqr+cptcovage;
+ /* for (k=1; k<=cptcovn;k++){ /\* Simple and product covariates without age* products *\/ */
+ for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */
+ cov[++ioffset]=covar[Tvar[k]][i];
+ }
+ for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitatives and Fixed covariates */
+ cov[++ioffset]=coqvar[iqv][i];
+ }
- /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]
- is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]
- has been calculated etc */
- /* For an individual i, wav[i] gives the number of effective waves */
- /* We compute the contribution to Likelihood of each effective transition
- mw[mi][i] is real wave of the mi th effectve wave */
- /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
- s2=s[mw[mi+1][i]][i];
- And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]
- But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
- meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
- */
- for(mi=1; mi<= wav[i]-1; mi++){
- for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */
- cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i];
- }
- for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */
- if(cotqvar[mw[mi][i]][iqtv][i] == -1){
- printf("i=%d, mi=%d, iqtv=%d, cotqvar[mw[mi][i]][iqtv][i]=%f",i,mi,iqtv,cotqvar[mw[mi][i]][iqtv][i]);
- }
- cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i];
- }
- /* ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; */
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact; /* Should be changed here */
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
- }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]
+ is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]
+ has been calculated etc */
+ /* For an individual i, wav[i] gives the number of effective waves */
+ /* We compute the contribution to Likelihood of each effective transition
+ mw[mi][i] is real wave of the mi th effectve wave */
+ /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
+ s2=s[mw[mi+1][i]][i];
+ And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]
+ But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
+ meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
+ */
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */
+ cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i];
+ }
+ for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */
+ if(cotqvar[mw[mi][i]][iqtv][i] == -1){
+ printf("i=%d, mi=%d, iqtv=%d, cotqvar[mw[mi][i]][iqtv][i]=%f",i,mi,iqtv,cotqvar[mw[mi][i]][iqtv][i]);
+ }
+ cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i];
+ }
+ /* ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; */
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact; /* Should be changed here */
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
- /* But now since version 0.9 we anticipate for bias at large stepm.
- * If stepm is larger than one month (smallest stepm) and if the exact delay
- * (in months) between two waves is not a multiple of stepm, we rounded to
- * the nearest (and in case of equal distance, to the lowest) interval but now
- * we keep into memory the bias bh[mi][i] and also the previous matrix product
- * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
- * probability in order to take into account the bias as a fraction of the way
- * from savm to out if bh is negative or even beyond if bh is positive. bh varies
- * -stepm/2 to stepm/2 .
- * For stepm=1 the results are the same as for previous versions of Imach.
- * For stepm > 1 the results are less biased than in previous versions.
- */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- bbh=(double)bh[mi][i]/(double)stepm;
- /* bias bh is positive if real duration
- * is higher than the multiple of stepm and negative otherwise.
- */
- /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
- if( s2 > nlstate){
- /* i.e. if s2 is a death state and if the date of death is known
- then the contribution to the likelihood is the probability to
- die between last step unit time and current step unit time,
- which is also equal to probability to die before dh
- minus probability to die before dh-stepm .
- In version up to 0.92 likelihood was computed
- as if date of death was unknown. Death was treated as any other
- health state: the date of the interview describes the actual state
- and not the date of a change in health state. The former idea was
- to consider that at each interview the state was recorded
- (healthy, disable or death) and IMaCh was corrected; but when we
- introduced the exact date of death then we should have modified
- the contribution of an exact death to the likelihood. This new
- contribution is smaller and very dependent of the step unit
- stepm. It is no more the probability to die between last interview
- and month of death but the probability to survive from last
- interview up to one month before death multiplied by the
- probability to die within a month. Thanks to Chris
- Jackson for correcting this bug. Former versions increased
- mortality artificially. The bad side is that we add another loop
- which slows down the processing. The difference can be up to 10%
- lower mortality.
- */
- /* If, at the beginning of the maximization mostly, the
- cumulative probability or probability to be dead is
- constant (ie = 1) over time d, the difference is equal to
- 0. out[s1][3] = savm[s1][3]: probability, being at state
- s1 at precedent wave, to be dead a month before current
- wave is equal to probability, being at state s1 at
- precedent wave, to be dead at mont of the current
- wave. Then the observed probability (that this person died)
- is null according to current estimated parameter. In fact,
- it should be very low but not zero otherwise the log go to
- infinity.
- */
+ /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+ /* But now since version 0.9 we anticipate for bias at large stepm.
+ * If stepm is larger than one month (smallest stepm) and if the exact delay
+ * (in months) between two waves is not a multiple of stepm, we rounded to
+ * the nearest (and in case of equal distance, to the lowest) interval but now
+ * we keep into memory the bias bh[mi][i] and also the previous matrix product
+ * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
+ * probability in order to take into account the bias as a fraction of the way
+ * from savm to out if bh is negative or even beyond if bh is positive. bh varies
+ * -stepm/2 to stepm/2 .
+ * For stepm=1 the results are the same as for previous versions of Imach.
+ * For stepm > 1 the results are less biased than in previous versions.
+ */
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ /* bias bh is positive if real duration
+ * is higher than the multiple of stepm and negative otherwise.
+ */
+ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
+ if( s2 > nlstate){
+ /* i.e. if s2 is a death state and if the date of death is known
+ then the contribution to the likelihood is the probability to
+ die between last step unit time and current step unit time,
+ which is also equal to probability to die before dh
+ minus probability to die before dh-stepm .
+ In version up to 0.92 likelihood was computed
+ as if date of death was unknown. Death was treated as any other
+ health state: the date of the interview describes the actual state
+ and not the date of a change in health state. The former idea was
+ to consider that at each interview the state was recorded
+ (healthy, disable or death) and IMaCh was corrected; but when we
+ introduced the exact date of death then we should have modified
+ the contribution of an exact death to the likelihood. This new
+ contribution is smaller and very dependent of the step unit
+ stepm. It is no more the probability to die between last interview
+ and month of death but the probability to survive from last
+ interview up to one month before death multiplied by the
+ probability to die within a month. Thanks to Chris
+ Jackson for correcting this bug. Former versions increased
+ mortality artificially. The bad side is that we add another loop
+ which slows down the processing. The difference can be up to 10%
+ lower mortality.
+ */
+ /* If, at the beginning of the maximization mostly, the
+ cumulative probability or probability to be dead is
+ constant (ie = 1) over time d, the difference is equal to
+ 0. out[s1][3] = savm[s1][3]: probability, being at state
+ s1 at precedent wave, to be dead a month before current
+ wave is equal to probability, being at state s1 at
+ precedent wave, to be dead at mont of the current
+ wave. Then the observed probability (that this person died)
+ is null according to current estimated parameter. In fact,
+ it should be very low but not zero otherwise the log go to
+ infinity.
+ */
/* #ifdef INFINITYORIGINAL */
/* lli=log(out[s1][s2] - savm[s1][s2]); */
/* #else */
/* else */
/* lli=log(out[s1][s2] - savm[s1][s2]); */
/* #endif */
- lli=log(out[s1][s2] - savm[s1][s2]);
+ lli=log(out[s1][s2] - savm[s1][s2]);
- } else if ( s2==-1 ) { /* alive */
- for (j=1,survp=0. ; j<=nlstate; j++)
- survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
- /*survp += out[s1][j]; */
- lli= log(survp);
- }
- else if (s2==-4) {
- for (j=3,survp=0. ; j<=nlstate; j++)
- survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
- lli= log(survp);
- }
- else if (s2==-5) {
- for (j=1,survp=0. ; j<=2; j++)
- survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
- lli= log(survp);
- }
- else{
- lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
- /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
- }
- /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
- /*if(lli ==000.0)*/
- /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- /* if (lli < log(mytinydouble)){ */
- /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */
- /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */
- /* } */
- } /* end of wave */
- } /* end of individual */
- } else if(mle==2){
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<=dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ } else if ( s2==-1 ) { /* alive */
+ for (j=1,survp=0. ; j<=nlstate; j++)
+ survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ /*survp += out[s1][j]; */
+ lli= log(survp);
+ }
+ else if (s2==-4) {
+ for (j=3,survp=0. ; j<=nlstate; j++)
+ survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ lli= log(survp);
+ }
+ else if (s2==-5) {
+ for (j=1,survp=0. ; j<=2; j++)
+ survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ lli= log(survp);
+ }
+ else{
+ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+ /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
+ }
+ /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
+ /*if(lli ==000.0)*/
+ /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ /* if (lli < log(mytinydouble)){ */
+ /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */
+ /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */
+ /* } */
+ } /* end of wave */
+ } /* end of individual */
+ } else if(mle==2){
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<=dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- bbh=(double)bh[mi][i]/(double)stepm;
- lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- } /* end of wave */
- } /* end of individual */
- } else if(mle==3){ /* exponential inter-extrapolation */
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ } /* end of wave */
+ } /* end of individual */
+ } else if(mle==3){ /* exponential inter-extrapolation */
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- bbh=(double)bh[mi][i]/(double)stepm;
- lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- } /* end of wave */
- } /* end of individual */
- }else if (mle==4){ /* ml=4 no inter-extrapolation */
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ } /* end of wave */
+ } /* end of individual */
+ }else if (mle==4){ /* ml=4 no inter-extrapolation */
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- if( s2 > nlstate){
- lli=log(out[s1][s2] - savm[s1][s2]);
- } else if ( s2==-1 ) { /* alive */
- for (j=1,survp=0. ; j<=nlstate; j++)
- survp += out[s1][j];
- lli= log(survp);
- }else{
- lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
- }
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ if( s2 > nlstate){
+ lli=log(out[s1][s2] - savm[s1][s2]);
+ } else if ( s2==-1 ) { /* alive */
+ for (j=1,survp=0. ; j<=nlstate; j++)
+ survp += out[s1][j];
+ lli= log(survp);
+ }else{
+ lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
+ }
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
/* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
- } /* end of wave */
- } /* end of individual */
- }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */
- for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
- for(mi=1; mi<= wav[i]-1; mi++){
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- }
+ } /* end of wave */
+ } /* end of individual */
+ }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */
+ for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ for(mi=1; mi<= wav[i]-1; mi++){
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+ }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
- ipmx +=1;
- sw += weight[i];
- ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/
- } /* end of wave */
- } /* end of individual */
- } /* End of if */
- for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
- /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
- l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
- return -l;
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
+ ipmx +=1;
+ sw += weight[i];
+ ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/
+ } /* end of wave */
+ } /* end of individual */
+ } /* End of if */
+ for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
+ /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
+ l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+ return -l;
}
/*************** log-likelihood *************/
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
ioffset=2+nagesqr+cptcovage;
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
- for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */
+ for (k=1; k<=ncoveff+nqfveff;k++){ /* Simple and product fixed covariates without age* products */
cov[++ioffset]=covar[Tvar[k]][i];
}
- for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitatives Fixed covariates */
- cov[++ioffset]=coqvar[iqv][i];
+ for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitative fixed covariates */
+ cov[++ioffset]=coqvar[Tvar[iqv]][i];
}
- for(mi=1; mi<= wav[i]-1; mi++){
+ for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */
for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */
cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i];
}
}
/************ Frequencies ********************/
- void freqsummary(char fileres[], 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 */
+void freqsummary(char fileres[], 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 */
- int i, m, jk, j1, bool, z1,j;
- int iind=0, iage=0;
- int mi; /* Effective wave */
- int first;
- double ***freq; /* Frequencies */
- double *meanq;
- double **meanqt;
- double *pp, **prop, *posprop, *pospropt;
- double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0;
- char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];
- double agebegin, ageend;
+ int i, m, jk, j1, bool, z1,j;
+ int iind=0, iage=0;
+ int mi; /* Effective wave */
+ int first;
+ double ***freq; /* Frequencies */
+ double *meanq;
+ double **meanqt;
+ double *pp, **prop, *posprop, *pospropt;
+ double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0;
+ char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];
+ double agebegin, ageend;
- pp=vector(1,nlstate);
- prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+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); */
- meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */
- meanqt=matrix(1,lastpass,1,nqtveff);
- strcpy(fileresp,"P_");
- strcat(fileresp,fileresu);
- /*strcat(fileresphtm,fileresu);*/
- if((ficresp=fopen(fileresp,"w"))==NULL) {
- printf("Problem with prevalence resultfile: %s\n", fileresp);
- fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
- exit(0);
- }
+ pp=vector(1,nlstate);
+ prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+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); */
+ meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */
+ meanqt=matrix(1,lastpass,1,nqtveff);
+ strcpy(fileresp,"P_");
+ strcat(fileresp,fileresu);
+ /*strcat(fileresphtm,fileresu);*/
+ if((ficresp=fopen(fileresp,"w"))==NULL) {
+ printf("Problem with prevalence resultfile: %s\n", fileresp);
+ fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
+ exit(0);
+ }
- strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
- if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
- printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
- fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
- fflush(ficlog);
- exit(70);
- }
- else{
- fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
+ strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
+ if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
+ printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
+ fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
+ fflush(ficlog);
+ exit(70);
+ }
+ else{
+ fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
- fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
- }
-
fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
+ fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+ }
+ fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
- strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
- if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
- printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
- fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
- fflush(ficlog);
- exit(70);
- }
- else{
- fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
+ strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
+ if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
+ printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
+ fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
+ fflush(ficlog);
+ exit(70);
+ }
+ else{
+ fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
- fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
- }
-
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
+ fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+ }
+ fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
- freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
- j1=0;
+ freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
+ j1=0;
- j=ncoveff;
- if (cptcovn<1) {j=1;ncodemax[1]=1;}
+ j=ncoveff;
+ if (cptcovn<1) {j=1;ncodemax[1]=1;}
- first=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
- med_educ V1=1 V2=0,
- 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
- */
+ /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
+ reference=low_education V1=0,V2=0
+ med_educ V1=1 V2=0,
+ 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 excluding varying and quantitatives */
- 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 (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination excluding varying and quantitatives */
+ 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++)
- prop[i][m]=0;
- posprop[i]=0;
- pospropt[i]=0;
- }
- for (z1=1; z1<= nqfveff; z1++) {
- meanq[z1]+=0.;
- for(m=1;m<=lastpass;m++){
- meanqt[m][z1]=0.;
- }
- }
+ for (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 comination of covariate j1, we count and print the frequencies */
- for (iind=1; iind<=imx; iind++) { /* For each individual iind */
- bool=1;
- if (nqfveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
- for (z1=1; z1<= nqfveff; z1++) {
- meanq[z1]+=coqvar[Tvar[z1]][iind];
- }
- for (z1=1; z1<=ncoveff; z1++) {
- /* if(Tvaraff[z1] ==-20){ */
- /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
- /* }else if(Tvaraff[z1] ==-10){ */
- /* /\* sumnew+=coqvar[z1][iind]; *\/ */
- /* }else */
- if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
- /* Tests if this individual i responded to j1 (V4=1 V3=0) */
- bool=0;
- /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",
- bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
- j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
- /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
- }
- } /* end z1 */
- } /* cptcovn > 0 */
-
- if (bool==1){ /* We selected an individual iin satisfying combination j1 */
- /* for(m=firstpass; m<=lastpass; m++){ */
- for(mi=1; mi<wav[iind];mi++){
- m=mw[mi][iind];
- /* 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 (m<lastpass) {
- /* if(s[m][iind]==4 && s[m+1][iind]==4) */
- /* printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
- if(s[m][iind]==-1)
- printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
- freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
- /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
- freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
- }
- }
- if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
- dateintsum=dateintsum+k2;
- k2cpt++;
- /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
- }
- /*}*/
- } /* end m */
- } /* end bool */
- } /* end iind = 1 to imx */
- /* prop[s][age] is feeded for any initial and valid live state as well as
- freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
-
-
- /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
- pstamp(ficresp);
- if (ncoveff>0) {
- fprintf(ficresp, "\n#********** Variable ");
- fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
- fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
- for (z1=1; z1<=ncoveff; z1++){
- fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- }
- fprintf(ficresp, "**********\n#");
- fprintf(ficresphtm, "**********</h3>\n");
- fprintf(ficresphtmfr, "**********</h3>\n");
- fprintf(ficlog, "\n#********** Variable ");
- for (z1=1; z1<=ncoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficlog, "**********\n");
- }
-
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
- for(i=1; i<=nlstate;i++) {
- fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
- fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
- }
- fprintf(ficresp, "\n");
- fprintf(ficresphtm, "\n");
+ dateintsum=0;
+ k2cpt=0;
+ /* For that comination of covariate j1, we count and print the frequencies */
+ for (iind=1; iind<=imx; iind++) { /* For each individual iind */
+ bool=1;
+ if (nqfveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+ for (z1=1; z1<= nqfveff; z1++) {
+ meanq[z1]+=coqvar[Tvar[z1]][iind];
+ }
+ for (z1=1; z1<=ncoveff; z1++) {
+ /* if(Tvaraff[z1] ==-20){ */
+ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
+ /* }else if(Tvaraff[z1] ==-10){ */
+ /* /\* sumnew+=coqvar[z1][iind]; *\/ */
+ /* }else */
+ if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
+ /* Tests if this individual i responded to j1 (V4=1 V3=0) */
+ bool=0;
+ /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",
+ bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
+ j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+ /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
+ }
+ } /* end z1 */
+ } /* cptcovn > 0 */
+
+ if (bool==1){ /* We selected an individual iin satisfying combination j1 */
+ /* for(m=firstpass; m<=lastpass; m++){ */
+ for(mi=1; mi<wav[iind];mi++){
+ m=mw[mi][iind];
+ /* 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 (m<lastpass) {
+ /* if(s[m][iind]==4 && s[m+1][iind]==4) */
+ /* printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
+ if(s[m][iind]==-1)
+ printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
+ freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
+ /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
+ freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
+ }
+ }
+ if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
+ dateintsum=dateintsum+k2;
+ k2cpt++;
+ /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
+ }
+ /*}*/
+ } /* end m */
+ } /* end bool */
+ } /* end iind = 1 to imx */
+ /* prop[s][age] is feeded for any initial and valid live state as well as
+ freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
+
+
+ /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
+ pstamp(ficresp);
+ if (ncoveff>0) {
+ fprintf(ficresp, "\n#********** Variable ");
+ fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
+ fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
+ for (z1=1; z1<=ncoveff; z1++){
+ fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ }
+ fprintf(ficresp, "**********\n#");
+ fprintf(ficresphtm, "**********</h3>\n");
+ fprintf(ficresphtmfr, "**********</h3>\n");
+ fprintf(ficlog, "\n#********** Variable ");
+ for (z1=1; z1<=ncoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficlog, "**********\n");
+ }
+ fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
+ for(i=1; i<=nlstate;i++) {
+ fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
+ fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
+ }
+ fprintf(ficresp, "\n");
+ fprintf(ficresphtm, "\n");
- /* Header of frequency table by age */
-
fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
- fprintf(ficresphtmfr,"<th>Age</th> ");
- for(jk=-1; jk <=nlstate+ndeath; jk++){
- for(m=-1; m <=nlstate+ndeath; m++){
- if(jk!=0 && m!=0)
- fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
- }
- }
- fprintf(ficresphtmfr, "\n");
+ /* Header of frequency table by age */
+ fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
+ fprintf(ficresphtmfr,"<th>Age</th> ");
+ for(jk=-1; jk <=nlstate+ndeath; jk++){
+ for(m=-1; m <=nlstate+ndeath; m++){
+ if(jk!=0 && m!=0)
+ fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
+ }
+ }
+ fprintf(ficresphtmfr, "\n");
- /* For each age */
- for(iage=iagemin; iage <= iagemax+3; iage++){
- fprintf(ficresphtm,"<tr>");
- if(iage==iagemax+1){
- fprintf(ficlog,"1");
- fprintf(ficresphtmfr,"<tr><th>0</th> ");
- }else if(iage==iagemax+2){
- fprintf(ficlog,"0");
- fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
- }else if(iage==iagemax+3){
- fprintf(ficlog,"Total");
- fprintf(ficresphtmfr,"<tr><th>Total</th> ");
- }else{
- if(first==1){
- first=0;
- printf("See log file for details...\n");
- }
- fprintf(ficresphtmfr,"<tr><th>%d</th> ",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);
- }
- }
-
- 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] */
- }
- for(jk=1; jk <=nlstate ; jk++){
- if(pos>=1.e-5){
- if(first==1)
- printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
- fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
- }else{
- if(first==1)
- printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
- fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
- }
- if( iage <= iagemax){
- if(pos>=1.e-5){
- fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
- fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",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,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
- }
- }
- pospropt[jk] +=posprop[jk];
- } /* end loop jk */
- /* pospropt=0.; */
- for(jk=-1; jk <=nlstate+ndeath; jk++){
- for(m=-1; m <=nlstate+ndeath; m++){
- if(freq[jk][m][iage] !=0 ) { /* minimizing output */
- if(first==1){
- printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
- }
- fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
- }
- if(jk!=0 && m!=0)
- fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
- }
- } /* end loop jk */
- posproptt=0.;
- for(jk=1; jk <=nlstate; jk++){
- posproptt += pospropt[jk];
- }
- fprintf(ficresphtmfr,"</tr>\n ");
- if(iage <= iagemax){
- fprintf(ficresp,"\n");
- fprintf(ficresphtm,"</tr>\n");
- }
- if(first==1)
- printf("Others in log...\n");
- fprintf(ficlog,"\n");
- } /* end loop age iage */
- fprintf(ficresphtm,"<tr><th>Tot</th>");
- for(jk=1; jk <=nlstate ; jk++){
- if(posproptt < 1.e-5){
- fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);
- }else{
- fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);
- }
- }
- fprintf(ficresphtm,"</tr>\n");
- fprintf(ficresphtm,"</table>\n");
- fprintf(ficresphtmfr,"</table>\n");
- if(posproptt < 1.e-5){
-
fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
-
fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",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 <p> This combination (%d) is valid and result will be produced.</p>",j1);
- invalidvarcomb[j1]=0;
- }
- fprintf(ficresphtmfr,"</table>\n");
- } /* end selected combination of covariate j1 */
- dateintmean=dateintsum/k2cpt;
+ /* For each age */
+ for(iage=iagemin; iage <= iagemax+3; iage++){
+ fprintf(ficresphtm,"<tr>");
+ if(iage==iagemax+1){
+ fprintf(ficlog,"1");
+ fprintf(ficresphtmfr,"<tr><th>0</th> ");
+ }else if(iage==iagemax+2){
+ fprintf(ficlog,"0");
+ fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
+ }else if(iage==iagemax+3){
+ fprintf(ficlog,"Total");
+ fprintf(ficresphtmfr,"<tr><th>Total</th> ");
+ }else{
+ if(first==1){
+ first=0;
+ printf("See log file for details...\n");
+ }
+ fprintf(ficresphtmfr,"<tr><th>%d</th> ",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);
+ }
+ }
+
+ 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] */
+ }
+ for(jk=1; jk <=nlstate ; jk++){
+ if(pos>=1.e-5){
+ if(first==1)
+ printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+ fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+ }else{
+ if(first==1)
+ printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+ fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+ }
+ if( iage <= iagemax){
+ if(pos>=1.e-5){
+ fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+ fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",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,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
+ }
+ }
+ pospropt[jk] +=posprop[jk];
+ } /* end loop jk */
+ /* pospropt=0.; */
+ for(jk=-1; jk <=nlstate+ndeath; jk++){
+ for(m=-1; m <=nlstate+ndeath; m++){
+ if(freq[jk][m][iage] !=0 ) { /* minimizing output */
+ if(first==1){
+ printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
+ }
+ fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
+ }
+ if(jk!=0 && m!=0)
+ fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
+ }
+ } /* end loop jk */
+ posproptt=0.;
+ for(jk=1; jk <=nlstate; jk++){
+ posproptt += pospropt[jk];
+ }
+ fprintf(ficresphtmfr,"</tr>\n ");
+ if(iage <= iagemax){
+ fprintf(ficresp,"\n");
+ fprintf(ficresphtm,"</tr>\n");
+ }
+ if(first==1)
+ printf("Others in log...\n");
+ fprintf(ficlog,"\n");
+ } /* end loop age iage */
+ fprintf(ficresphtm,"<tr><th>Tot</th>");
+ for(jk=1; jk <=nlstate ; jk++){
+ if(posproptt < 1.e-5){
+ fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);
+ }else{
+ fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);
+ }
+ }
+ fprintf(ficresphtm,"</tr>\n");
+ fprintf(ficresphtm,"</table>\n");
+ fprintf(ficresphtmfr,"</table>\n");
+ if(posproptt < 1.e-5){
+ fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
+ fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",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 <p> This combination (%d) is valid and result will be produced.</p>",j1);
+ invalidvarcomb[j1]=0;
+ }
+ fprintf(ficresphtmfr,"</table>\n");
+ } /* end selected combination of covariate j1 */
+ dateintmean=dateintsum/k2cpt;
- fclose(ficresp);
- fclose(ficresphtm);
- 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(pospropt,1,nlstate);
- free_vector(posprop,1,nlstate);
- free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);
- free_vector(pp,1,nlstate);
- /* End of freqsummary */
- }
+ fclose(ficresp);
+ fclose(ficresphtm);
+ 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(pospropt,1,nlstate);
+ free_vector(posprop,1,nlstate);
+ free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);
+ free_vector(pp,1,nlstate);
+ /* End of freqsummary */
+}
/************ Prevalence ********************/
void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
return 1;
}
coqvar[iv][i]=dval;
+ covar[ncovcol+iv][i]=dval; /* including qvar in standard covar for performance reasons */
}
strcpy(line,stra);
}/* end loop nqv */
cptcovprod--;
cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
- Typevar[k]=1; /* 2 for age product */
+ Typevar[k]=1; /* 1 for age product */
cptcovage++; /* Sums the number of covariates which include age as a product */
Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
/*printf("stre=%s ", stre);*/
/* Decodemodel knows only the grammar (simple, product, age*) of the model but not what kind
of variable (dummy vs quantitative, fixed vs time varying) is behind */
-/* ncovcol= 1, nqv=1, ntv=2, nqtv= 1 = 5 possible variables data
- model= V2 + V4 +V3 + V4*V3 + V5*age + V5 , V1 is not used saving its place
- k = 1 2 3 4 5 6
- Tvar[k]= 2 4 3 1+1+2+1+1=6 5 5
- Typevar[k]=0 0 0 2 1 0
+/* ncovcol= 1, nqv=1 | ntv=2, nqtv= 1 = 5 possible variables data: 2 fixed 3, varying
+ model= V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place
+ k = 1 2 3 4 5 6 7 8 9
+ Tvar[k]= 5 4 3 1+1+2+1+1=6 5 2 7 1 5
+ Typevar[k]= 0 0 0 2 1 0 2 1 1
+ Fixed[Tvar[k]]1 1 1 1 2 0 1 2 3
+ Dummy[Tvar[k]]1 0 0 0 2 1 1 2 3
*/
/* Dispatching between quantitative and time varying covariates */
+ /* If Tvar[k] >ncovcol it is a product */
/* Tvar[k] is the value n of Vn with n varying for 1 to nvcol, or p Vp=Vn*Vm for product */
+ /* Computing effective variables, ie used by the model, that is from the cptcovt variables */
for(k=1, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
- if (Tvar[k] <=ncovcol){ /* Simple fixed dummy covariatee */
+ if (Tvar[k] <=ncovcol && (Typevar[k]==0 || Typevar[k]==2)){ /* Simple or product fixed dummy covariatee */
+ Fixed[Tvar[k]]= 0;
+ Dummy[Tvar[k]]= 0;
ncoveff++;
}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/
+ Fixed[Tvar[k]]= 0;
+ Dummy[Tvar[k]]= 1;
nqfveff++; /* Only simple fixed quantitative variable */
}else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 0;
ntveff++; /* Only simple time varying dummy variable */
- }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){
- nqtveff++;/* Only simple time varying quantitative variable */
+ }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
+ if( Typevar[k]==0){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ nqtveff++;/* Only simple time varying quantitative variable */
+ }
+ }else if (Typevar[k] == 2) {
+ for(k1=1; k1 <= cptcovprodnoage; k1++){
+ if(Tvard[k1][1] <=ncovcol){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 0;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[Tvar[k]]= 0;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 0;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }
+ }else if(Tvard[k1][1] <=ncovcol+nqv){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[Tvar[k]]= 0;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[Tvar[k]]= 0;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }
+ }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 0;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }
+ }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[Tvar[k]]= 1;
+ Dummy[Tvar[k]]= 1;
+ }
+ }else{
+ printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+ fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+ }
+ } /* end k1 */
}else{
- printf("Other types in effective covariates \n");
+ printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
+ fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
}
+ printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[Tvar[k]],Dummy[Tvar[k]]);
+ fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[Tvar[k]],Dummy[Tvar[k]]);
}
-
+ printf("Model=%s\n\
+Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\
+Fixed[Tvar[k]] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product \n\
+Dummy[Tvar[k]] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
+
printf("ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn);
fprintf(ficlog,"ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn);
return (0); /* with covar[new additional covariate if product] and Tage if age */
/* Scans npar lines */
for(i=1; i <=npar; i++){
- count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk);
+ count=fscanf(ficpar,"%1d%1d%d",&i1,&j1,&jk);
if(count != 3){
- printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
+ printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
- fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
+ fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
- exit(1);
+ exit(1);
}else{
- if(mle==1)
- printf("%1d%1d%1d",i1,j1,jk);
- }
- fprintf(ficlog,"%1d%1d%1d",i1,j1,jk);
- fprintf(ficparo,"%1d%1d%1d",i1,j1,jk);
+ if(mle==1)
+ printf("%1d%1d%d",i1,j1,jk);
+ }
+ fprintf(ficlog,"%1d%1d%d",i1,j1,jk);
+ fprintf(ficparo,"%1d%1d%d",i1,j1,jk);
for(j=1; j <=i; j++){
- fscanf(ficpar," %le",&matcov[i][j]);
- if(mle==1){
- printf(" %.5le",matcov[i][j]);
- }
- fprintf(ficlog," %.5le",matcov[i][j]);
- fprintf(ficparo," %.5le",matcov[i][j]);
+ fscanf(ficpar," %le",&matcov[i][j]);
+ if(mle==1){
+ printf(" %.5le",matcov[i][j]);
+ }
+ fprintf(ficlog," %.5le",matcov[i][j]);
+ fprintf(ficparo," %.5le",matcov[i][j]);
}
fscanf(ficpar,"\n");
numlinepar++;
/* End of read covariance matrix npar lines */
for(i=1; i <=npar; i++)
for(j=i+1;j<=npar;j++)
- matcov[i][j]=matcov[j][i];
+ matcov[i][j]=matcov[j][i];
if(mle==1)
printf("\n");
k=1 Tvar[1]=2 (from V2)
*/
Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
- Typevar=ivector(1,NCOVMAX); /* -1 to 4 */
+ Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */
+ Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */
+ Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */
/* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs).
For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4,
Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
for (i=1; i<=imx; i++){
dcwave[i]=-1;
for (m=firstpass; m<=lastpass; m++)
- if (s[m][i]>nlstate) {
- dcwave[i]=m;
- /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
- break;
- }
+ if (s[m][i]>nlstate) {
+ dcwave[i]=m;
+ /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
+ break;
+ }
}
-
+
for (i=1; i<=imx; i++) {
if (wav[i]>0){
- ageexmed[i]=agev[mw[1][i]][i];
- j=wav[i];
- agecens[i]=1.;
-
- if (ageexmed[i]> 1 && wav[i] > 0){
- agecens[i]=agev[mw[j][i]][i];
- cens[i]= 1;
- }else if (ageexmed[i]< 1)
- cens[i]= -1;
- if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
- cens[i]=0 ;
+ ageexmed[i]=agev[mw[1][i]][i];
+ j=wav[i];
+ agecens[i]=1.;
+
+ if (ageexmed[i]> 1 && wav[i] > 0){
+ agecens[i]=agev[mw[j][i]][i];
+ cens[i]= 1;
+ }else if (ageexmed[i]< 1)
+ cens[i]= -1;
+ if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
+ cens[i]=0 ;
}
else cens[i]=-1;
}
for (i=1;i<=NDIM;i++) {
for (j=1;j<=NDIM;j++)
- ximort[i][j]=(i == j ? 1.0 : 0.0);
+ ximort[i][j]=(i == j ? 1.0 : 0.0);
}
/*p[1]=0.0268; p[NDIM]=0.083;*/
free_ivector(ncodemax,1,NCOVMAX);
free_ivector(ncodemaxwundef,1,NCOVMAX);
+ free_ivector(Dummy,-1,NCOVMAX);
+ free_ivector(Fixed,-1,NCOVMAX);
free_ivector(Typevar,-1,NCOVMAX);
free_ivector(Tvar,1,NCOVMAX);
free_ivector(Tprod,1,NCOVMAX);
git://henry.ined.fr / .git / commitdiff