/* $Id$
$State$
$Log$
+ Revision 1.226 2016/07/12 18:42:34 brouard
+ Summary: temp
+
Revision 1.225 2016/07/12 08:40:03 brouard
Summary: saving but not running
Short summary of the programme:
- This program computes Healthy Life Expectancies from
- cross-longitudinal data. Cross-longitudinal data consist in: -1- a
- first survey ("cross") where individuals from different ages are
- interviewed on their health status or degree of disability (in the
- case of a health survey which is our main interest) -2- at least a
- second wave of interviews ("longitudinal") which measure each change
- (if any) in individual health status. Health expectancies are
- computed from the time spent in each health state according to a
- model. More health states you consider, more time is necessary to reach the
- Maximum Likelihood of the parameters involved in the model. The
- simplest model is the multinomial logistic model where pij is the
- probability to be observed in state j at the second wave
- conditional to be observed in state i at the first wave. Therefore
- the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
- 'age' is age and 'sex' is a covariate. If you want to have a more
- complex model than "constant and age", you should modify the program
- where the markup *Covariates have to be included here again* invites
- you to do it. More covariates you add, slower the
+ This program computes Healthy Life Expectancies or State-specific
+ (if states aren't health statuses) Expectancies from
+ cross-longitudinal data. Cross-longitudinal data consist in:
+
+ -1- a first survey ("cross") where individuals from different ages
+ are interviewed on their health status or degree of disability (in
+ the case of a health survey which is our main interest)
+
+ -2- at least a second wave of interviews ("longitudinal") which
+ measure each change (if any) in individual health status. Health
+ expectancies are computed from the time spent in each health state
+ according to a model. More health states you consider, more time is
+ necessary to reach the Maximum Likelihood of the parameters involved
+ in the model. The simplest model is the multinomial logistic model
+ where pij is the probability to be observed in state j at the second
+ wave conditional to be observed in state i at the first
+ wave. Therefore the model is: log(pij/pii)= aij + bij*age+ cij*sex +
+ etc , where 'age' is age and 'sex' is a covariate. If you want to
+ have a more complex model than "constant and age", you should modify
+ the program where the markup *Covariates have to be included here
+ again* invites you to do it. More covariates you add, slower the
convergence.
The advantage of this computer programme, compared to a simple
of the life expectancies. It also computes the period (stable) prevalence.
Back prevalence and projections:
- - back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj)
- Computes the back prevalence limit for any combination of covariate values k
- at any age between ageminpar and agemaxpar and returns it in **bprlim. In the loops,
- - **bprevalim(**bprlim, ***mobaverage, nlstate, *p, age, **oldm, **savm, **dnewm, **doldm, **dsavm, ftolpl, ncvyearp, k);
- - hBijx Back Probability to be in state i at age x-h being in j at x
+
+ - back_prevalence_limit(double *p, double **bprlim, double ageminpar,
+ double agemaxpar, double ftolpl, int *ncvyearp, double
+ dateprev1,double dateprev2, int firstpass, int lastpass, int
+ mobilavproj)
+
+ Computes the back prevalence limit for any combination of
+ covariate values k at any age between ageminpar and agemaxpar and
+ returns it in **bprlim. In the loops,
+
+ - **bprevalim(**bprlim, ***mobaverage, nlstate, *p, age, **oldm,
+ **savm, **dnewm, **doldm, **dsavm, ftolpl, ncvyearp, k);
+
+ - hBijx Back Probability to be in state i at age x-h being in j at x
Computes for any combination of covariates k and any age between bage and fage
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
- - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+
+ - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
Computes the transition matrix starting at age 'age' over
'nhstepm*hstepm*stepm' months (i.e. until
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
- nhstepm*hstepm matrices. Returns p3mat[i][j][h] after calling
- p3mat[i][j][h]=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\
- 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+ nhstepm*hstepm matrices.
+
+ Returns p3mat[i][j][h] after calling
+ p3mat[i][j][h]=matprod2(newm,
+ bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm,
+ dsavm,ij),\ 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
+ oldm);
Important routines
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */
-int *Fixed; /** Fixed[Tvar[k]] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
-int *Dummy; /** Dummy[Tvar[k]] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
+int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
+int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
int *Tage;
+int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */
+int *Tmodelind; /** Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/
int *Ndum; /** Freq of modality (tricode */
/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
-int **Tvard, *Tprod, cptcovprod, *Tvaraff, *invalidvarcomb;
+int **Tvard;
+int *Tprod;/**< Gives the k position of the k1 product */
+int *Tposprod; /**< Gives the k1 product from the k position */
+/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3
+ if V2+V1+V1*V4+age*V3+V3*V2 TProd[k1=2]=5 (V3*V2)
+ Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2
+*/
+int cptcovprod, *Tvaraff, *invalidvarcomb;
double *lsurv, *lpop, *tpop;
double ftol=FTOL; /**< Tolerance for computing Max Likelihood */
agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
cov[2]=agexact;
if(nagesqr==1)
- cov[3]= agexact*agexact;
+ cov[3]= agexact*agexact;
for (k=1; k<=cptcovn;k++)
- cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
- /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
- /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
- cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
- /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
-
-
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+
+
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
/* right multiplication of oldm by the current matrix */
int firstpass, int lastpass, int stepm, int weightopt, char model[])
{ /* Some frequencies */
- int i, m, jk, j1, bool, z1,j;
+ int i, m, jk, j1, bool, z1,j, k, iv;
int iind=0, iage=0;
int mi; /* Effective wave */
int first;
freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
j1=0;
- j=ncoveff;
+ /* j=ncoveff; /\* Only fixed dummy covariates *\/ */
+ j=cptcoveff; /* Only dummy covariates of the model */
if (cptcovn<1) {j=1;ncodemax[1]=1;}
first=1;
Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff
*/
- for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination excluding varying and quantitatives */
+ for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */
posproptt=0.;
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/
posprop[i]=0;
pospropt[i]=0;
}
- for (z1=1; z1<= nqfveff; z1++) {
- meanq[z1]+=0.;
- for(m=1;m<=lastpass;m++){
- meanqt[m][z1]=0.;
- }
- }
+ /* for (z1=1; z1<= nqfveff; z1++) { */
+ /* meanq[z1]+=0.; */
+ /* for(m=1;m<=lastpass;m++){ */
+ /* meanqt[m][z1]=0.; */
+ /* } */
+ /* } */
dateintsum=0;
k2cpt=0;
- /* For that comination of covariate j1, we count and print the frequencies */
+ /* For that combination of covariate j1, we count and print the frequencies in one pass */
for (iind=1; iind<=imx; iind++) { /* For each individual iind */
bool=1;
- if (nqfveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
- for (z1=1; z1<= nqfveff; z1++) {
- meanq[z1]+=coqvar[Tvar[z1]][iind];
- }
- for (z1=1; z1<=ncoveff; z1++) {
- /* if(Tvaraff[z1] ==-20){ */
- /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
- /* }else if(Tvaraff[z1] ==-10){ */
- /* /\* sumnew+=coqvar[z1][iind]; *\/ */
- /* }else */
- if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
- /* Tests if this individual i responded to j1 (V4=1 V3=0) */
- bool=0;
- /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",
- bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
- j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
- /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
- }
- } /* end z1 */
- } /* cptcovn > 0 */
-
- if (bool==1){ /* We selected an individual iin satisfying combination j1 */
+ if(anyvaryingduminmodel==0){ /* If All fixed covariates */
+ if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+ /* for (z1=1; z1<= nqfveff; z1++) { */
+ /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */
+ /* } */
+ for (z1=1; z1<=cptcoveff; z1++) {
+ /* if(Tvaraff[z1] ==-20){ */
+ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
+ /* }else if(Tvaraff[z1] ==-10){ */
+ /* /\* sumnew+=coqvar[z1][iind]; *\/ */
+ /* }else */
+ if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
+ /* Tests if this individual iind responded to j1 (V4=1 V3=0) */
+ bool=0;
+ /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",
+ bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
+ j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+ /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
+ } /* Onlyf fixed */
+ } /* end z1 */
+ } /* cptcovn > 0 */
+ } /* end any */
+ if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */
/* for(m=firstpass; m<=lastpass; m++){ */
- for(mi=1; mi<wav[iind];mi++){
+ for(mi=1; mi<wav[iind];mi++){ /* For that wave */
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(anyvaryingduminmodel==1){ /* Some are varying covariates */
+ for (z1=1; z1<=cptcoveff; z1++) {
+ if( Fixed[Tmodelind[z1]]==1){
+ iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+ if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */
+ bool=0;
+ }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
+ if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
+ bool=0;
+ }
+ }
}
- }
- 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); */
- }
- /*}*/
+ }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */
+ /* bool =0 we keep that guy which corresponds to the combination of dummy values */
+ if(bool==1){
+ /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
+ and mw[mi+1][iind]. dh depends on stepm. */
+ agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
+ ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
+ if(m >=firstpass && m <=lastpass){
+ k2=anint[m][iind]+(mint[m][iind]/12.);
+ /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
+ if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */
+ if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */
+ if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */
+ prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
+ if (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 */
+ }
+ } /* end if between passes */
+ 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 bool 2 */
} /* end m */
} /* end bool */
} /* end iind = 1 to imx */
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
pstamp(ficresp);
- if (ncoveff>0) {
+ /* if (ncoveff>0) { */
+ if (cptcoveff>0) {
fprintf(ficresp, "\n#********** Variable ");
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
- for (z1=1; z1<=ncoveff; z1++){
+ for (z1=1; z1<=cptcoveff; z1++){
fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(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)]);
+ for (z1=1; z1<=cptcoveff; 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\">");
}
/************ Prevalence ********************/
- void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
- {
- /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
- in each health status at the date of interview (if between dateprev1 and dateprev2).
- We still use firstpass and lastpass as another selection.
- */
+void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
+{
+ /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+ in each health status at the date of interview (if between dateprev1 and dateprev2).
+ We still use firstpass and lastpass as another selection.
+ */
- int i, m, jk, j1, bool, z1,j;
- int mi; /* Effective wave */
- int iage;
- double agebegin, ageend;
-
- double **prop;
- double posprop;
- double y2; /* in fractional years */
- int iagemin, iagemax;
- int first; /** to stop verbosity which is redirected to log file */
-
- iagemin= (int) agemin;
- iagemax= (int) agemax;
- /*pp=vector(1,nlstate);*/
- prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
- /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
- j1=0;
+ int i, m, jk, j1, bool, z1,j, iv;
+ int mi; /* Effective wave */
+ int iage;
+ double agebegin, ageend;
+
+ double **prop;
+ double posprop;
+ double y2; /* in fractional years */
+ int iagemin, iagemax;
+ int first; /** to stop verbosity which is redirected to log file */
+
+ iagemin= (int) agemin;
+ iagemax= (int) agemax;
+ /*pp=vector(1,nlstate);*/
+ prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
+ /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
+ j1=0;
- /*j=cptcoveff;*/
- if (cptcovn<1) {j=1;ncodemax[1]=1;}
+ /*j=cptcoveff;*/
+ if (cptcovn<1) {j=1;ncodemax[1]=1;}
- first=1;
- for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
- for (i=1; i<=nlstate; i++)
- for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
- prop[i][iage]=0.0;
+ first=1;
+ for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
+ for (i=1; i<=nlstate; i++)
+ for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
+ prop[i][iage]=0.0;
+ printf("Prevalence combination of varying and fixed dummies %d\n",j1);
+ /* fprintf(ficlog," V%d=%d ",Tvaraff[j1],nbcode[Tvaraff[j1]][codtabm(k,j1)]); */
+ fprintf(ficlog,"Prevalence combination of varying and fixed dummies %d\n",j1);
- for (i=1; i<=imx; i++) { /* Each individual */
- bool=1;
- if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
- for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
- if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
- bool=0;
- }
- if (bool==1) { /* For this combination of covariates values, this individual fits */
- /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
- for(mi=1; mi<wav[i];mi++){
- m=mw[mi][i];
- agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
- /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
- if(m >=firstpass && m <=lastpass){
- y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
- if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
- if(agev[m][i]==0) agev[m][i]=iagemax+1;
- if(agev[m][i]==1) agev[m][i]=iagemax+2;
- if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
- printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m);
- exit(1);
- }
- if (s[m][i]>0 && s[m][i]<=nlstate) {
- /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
- prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
- prop[s[m][i]][iagemax+3] += weight[i];
- } /* end valid statuses */
- } /* end selection of dates */
- } /* end selection of waves */
- } /* end effective waves */
- } /* end bool */
- }
- for(i=iagemin; i <= iagemax+3; i++){
- for(jk=1,posprop=0; jk <=nlstate ; jk++) {
- posprop += prop[jk][i];
- }
+ for (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<wav[i];mi++){ /* For this wave too look where individual can be counted V4=0 V3=0 */
+ m=mw[mi][i];
+ /* Tmodelind[z1]=k is the position of the varying covariate in the model, but which # within 1 to ntv? */
+ /* Tvar[Tmodelind[z1]] is the n of Vn; n-ncovcol-nqv is the first time varying covariate or iv */
+ for (z1=1; z1<=cptcoveff; z1++){
+ if( Fixed[Tmodelind[z1]]==1){
+ iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+ if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */
+ bool=0;
+ }else if( Fixed[Tmodelind[z1]]== 0) /* fixed */
+ if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
+ bool=0;
+ }
+ }
+ if(bool==1){ /* Otherwise we skip that wave/person */
+ agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
+ /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
+ if(m >=firstpass && m <=lastpass){
+ y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
+ if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
+ if(agev[m][i]==0) agev[m][i]=iagemax+1;
+ if(agev[m][i]==1) agev[m][i]=iagemax+2;
+ if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
+ printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m);
+ exit(1);
+ }
+ if (s[m][i]>0 && s[m][i]<=nlstate) {
+ /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
+ prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
+ prop[s[m][i]][iagemax+3] += weight[i];
+ } /* end valid statuses */
+ } /* end selection of dates */
+ } /* end selection of waves */
+ } /* end bool */
+ } /* end wave */
+ } /* end individual */
+ for(i=iagemin; i <= iagemax+3; i++){
+ for(jk=1,posprop=0; jk <=nlstate ; jk++) {
+ posprop += prop[jk][i];
+ }
- for(jk=1; jk <=nlstate ; jk++){
- if( i <= iagemax){
- if(posprop>=1.e-5){
- probs[i][jk][j1]= prop[jk][i]/posprop;
- } else{
- if(first==1){
- first=0;
- printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
- }
- }
- }
- }/* end jk */
- }/* end i */
+ for(jk=1; jk <=nlstate ; jk++){
+ if( i <= iagemax){
+ if(posprop>=1.e-5){
+ probs[i][jk][j1]= prop[jk][i]/posprop;
+ } else{
+ if(first==1){
+ first=0;
+ printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,j1,probs[i][jk][j1]);
+ }
+ }
+ }
+ }/* end jk */
+ }/* end i */
/*} *//* end i1 */
- } /* end j1 */
+ } /* end j1 */
- /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
- /*free_vector(pp,1,nlstate);*/
- free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE);
- } /* End of prevalence */
+ /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
+ /*free_vector(pp,1,nlstate);*/
+ free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE);
+} /* End of prevalence */
/************* Waves Concatenation ***************/
mw[mi][i] is the mi (mi=1 to wav[i]) effective wave of individual i
dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
and mw[mi+1][i]. dh depends on stepm.
- */
+ */
int i=0, mi=0, m=0, mli=0;
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
/* Treating live states */
for(i=1; i<=imx; i++){ /* For simple cases and if state is death */
mi=0; /* First valid wave */
- mli=0; /* Last valid wave */
+ mli=0; /* Last valid wave */
m=firstpass;
while(s[m][i] <= nlstate){ /* a live state */
- if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */
- mli=m-1;/* mw[++mi][i]=m-1; */
- }else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */
- mw[++mi][i]=m;
- mli=m;
+ if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */
+ mli=m-1;/* mw[++mi][i]=m-1; */
+ }else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */
+ mw[++mi][i]=m;
+ mli=m;
} /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */
if(m < lastpass){ /* m < lastpass, standard case */
- m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */
+ m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */
}
- else{ /* m >= lastpass, eventual special issue with warning */
+ else{ /* m >= lastpass, eventual special issue with warning */
#ifdef UNKNOWNSTATUSNOTCONTRIBUTING
- break;
+ break;
#else
- if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
- if(firsthree == 0){
- printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
- firsthree=1;
- }
- fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
- mw[++mi][i]=m;
- mli=m;
- }
- if(s[m][i]==-2){ /* Vital status is really unknown */
- nbwarn++;
- if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */
- printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
- fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
- }
- break;
- }
- break;
+ if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
+ if(firsthree == 0){
+ printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+ firsthree=1;
+ }
+ fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+ mw[++mi][i]=m;
+ mli=m;
+ }
+ if(s[m][i]==-2){ /* Vital status is really unknown */
+ nbwarn++;
+ if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */
+ printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+ fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+ }
+ break;
+ }
+ break;
#endif
- }/* End m >= lastpass */
+ }/* End m >= lastpass */
}/* end while */
- /* mi is the last effective wave, m is lastpass, mw[j][i] gives the # of j-th effective wave for individual i */
+ /* mi is the last effective wave, m is lastpass, mw[j][i] gives the # of j-th effective wave for individual i */
/* After last pass */
/* Treating death states */
if (s[m][i] > nlstate){ /* In a death state */
- /* if( mint[m][i]==mdc[m][i] && anint[m][i]==andc[m][i]){ /\* same date of death and date of interview *\/ */
- /* } */
+ /* if( mint[m][i]==mdc[m][i] && anint[m][i]==andc[m][i]){ /\* same date of death and date of interview *\/ */
+ /* } */
mi++; /* Death is another wave */
/* if(mi==0) never been interviewed correctly before death */
- /* Only death is a correct wave */
+ /* Only death is a correct wave */
mw[mi][i]=m;
}
#ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE
- else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */
+ else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */
/* m++; */
/* mi++; */
/* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */
/* mw[mi][i]=m; */
if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
- if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */
- nbwarn++;
- if(firstfiv==0){
- printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m );
- firstfiv=1;
- }else{
- fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m );
- }
- }else{ /* Death occured afer last wave potential bias */
- nberr++;
- if(firstwo==0){
- printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
- firstwo=1;
- }
- fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
- }
+ if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */
+ nbwarn++;
+ if(firstfiv==0){
+ printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m );
+ firstfiv=1;
+ }else{
+ fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m );
+ }
+ }else{ /* Death occured afer last wave potential bias */
+ nberr++;
+ if(firstwo==0){
+ printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ firstwo=1;
+ }
+ fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ }
}else{ /* end date of interview is known */
- /* death is known but not confirmed by death status at any wave */
- if(firstfour==0){
- printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
- firstfour=1;
- }
- fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ /* death is known but not confirmed by death status at any wave */
+ if(firstfour==0){
+ printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ firstfour=1;
+ }
+ fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
}
} /* end if date of death is known */
#endif
if(mi==0){
nbwarn++;
if(first==0){
- printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
- first=1;
+ printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
+ first=1;
}
if(first==1){
- fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
+ fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
}
} /* end mi==0 */
} /* End individuals */
for(i=1; i<=imx; i++){
for(mi=1; mi<wav[i];mi++){
if (stepm <=0)
- dh[mi][i]=1;
+ dh[mi][i]=1;
else{
- if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
- if (agedc[i] < 2*AGESUP) {
- j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);
- if(j==0) j=1; /* Survives at least one month after exam */
- else if(j<0){
- nberr++;
- printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- j=1; /* Temporary Dangerous patch */
- printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
- fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
- }
- k=k+1;
- if (j >= jmax){
- jmax=j;
- ijmax=i;
- }
- if (j <= jmin){
- jmin=j;
- ijmin=i;
- }
- sum=sum+j;
- /*if (j<0) printf("j=%d num=%d \n",j,i);*/
- /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
- }
- }
- else{
- j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
+ if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
+ if (agedc[i] < 2*AGESUP) {
+ j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);
+ if(j==0) j=1; /* Survives at least one month after exam */
+ else if(j<0){
+ nberr++;
+ printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ j=1; /* Temporary Dangerous patch */
+ printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
+ fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
+ }
+ k=k+1;
+ if (j >= jmax){
+ jmax=j;
+ ijmax=i;
+ }
+ if (j <= jmin){
+ jmin=j;
+ ijmin=i;
+ }
+ sum=sum+j;
+ /*if (j<0) printf("j=%d num=%d \n",j,i);*/
+ /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
+ }
+ }
+ else{
+ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
- k=k+1;
- if (j >= jmax) {
- jmax=j;
- ijmax=i;
- }
- else if (j <= jmin){
- jmin=j;
- ijmin=i;
- }
- /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
- /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
- if(j<0){
- nberr++;
- printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- }
- sum=sum+j;
- }
- jk= j/stepm;
- jl= j -jk*stepm;
- ju= j -(jk+1)*stepm;
- if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
- if(jl==0){
- dh[mi][i]=jk;
- bh[mi][i]=0;
- }else{ /* We want a negative bias in order to only have interpolation ie
- * to avoid the price of an extra matrix product in likelihood */
- dh[mi][i]=jk+1;
- bh[mi][i]=ju;
- }
- }else{
- if(jl <= -ju){
- dh[mi][i]=jk;
- bh[mi][i]=jl; /* bias is positive if real duration
- * is higher than the multiple of stepm and negative otherwise.
- */
- }
- else{
- dh[mi][i]=jk+1;
- bh[mi][i]=ju;
- }
- if(dh[mi][i]==0){
- dh[mi][i]=1; /* At least one step */
- bh[mi][i]=ju; /* At least one step */
- /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
- }
- } /* end if mle */
+ k=k+1;
+ if (j >= jmax) {
+ jmax=j;
+ ijmax=i;
+ }
+ else if (j <= jmin){
+ jmin=j;
+ ijmin=i;
+ }
+ /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
+ /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
+ if(j<0){
+ nberr++;
+ printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ }
+ sum=sum+j;
+ }
+ jk= j/stepm;
+ jl= j -jk*stepm;
+ ju= j -(jk+1)*stepm;
+ if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
+ if(jl==0){
+ dh[mi][i]=jk;
+ bh[mi][i]=0;
+ }else{ /* We want a negative bias in order to only have interpolation ie
+ * to avoid the price of an extra matrix product in likelihood */
+ dh[mi][i]=jk+1;
+ bh[mi][i]=ju;
+ }
+ }else{
+ if(jl <= -ju){
+ dh[mi][i]=jk;
+ bh[mi][i]=jl; /* bias is positive if real duration
+ * is higher than the multiple of stepm and negative otherwise.
+ */
+ }
+ else{
+ dh[mi][i]=jk+1;
+ bh[mi][i]=ju;
+ }
+ if(dh[mi][i]==0){
+ dh[mi][i]=1; /* At least one step */
+ bh[mi][i]=ju; /* At least one step */
+ /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
+ }
+ } /* end if mle */
}
} /* end wave */
}
jmean=sum/k;
printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);
fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %d) Max=%d (%d) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
- }
+}
/*********** Tricode ****************************/
void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum)
/* Loop on covariates without age and products and no quantitative variable */
/* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
- for (j=1; j<=cptcovsnq; j++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
- for (k=-1; k < maxncov; k++) Ndum[k]=0;
- for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the
- modality of this covariate Vj*/
- switch(Typevar[j]) {
- case 1: /* A real fixed dummy covariate */
- ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
- * If product of Vn*Vm, still boolean *:
- * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
- * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */
- /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
- modality of the nth covariate of individual i. */
- if (ij > modmaxcovj)
- modmaxcovj=ij;
- else if (ij < modmincovj)
- modmincovj=ij;
- if ((ij < -1) && (ij > NCOVMAX)){
- printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
- exit(1);
- }else
- Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
- /* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
- /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
- /* getting the maximum value of the modality of the covariate
- (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
- female ies 1, then modmaxcovj=1.*/
- break;
- case 2:
- break;
-
- }
- } /* end for loop on individuals i */
- printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
- fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
- cptcode=modmaxcovj;
- /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
- /*for (i=0; i<=cptcode; i++) {*/
- for (k=modmincovj; k<=modmaxcovj; k++) { /* k=-1 ? 0 and 1*//* For each value k of the modality of model-cov j */
- printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
- fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
- if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */
- if( k != -1){
- ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th
- covariate for which somebody answered excluding
- undefined. Usually 2: 0 and 1. */
+ for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
+ for (j=-1; (j < maxncov); j++) Ndum[j]=0;
+ if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */
+ switch(Fixed[k]) {
+ case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
+ for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/
+ ij=(int)(covar[Tvar[k]][i]);
+ /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
+ * If product of Vn*Vm, still boolean *:
+ * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
+ * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */
+ /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
+ modality of the nth covariate of individual i. */
+ if (ij > modmaxcovj)
+ modmaxcovj=ij;
+ else if (ij < modmincovj)
+ modmincovj=ij;
+ if ((ij < -1) && (ij > NCOVMAX)){
+ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
+ exit(1);
+ }else
+ Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
+ /* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
+ /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
+ /* getting the maximum value of the modality of the covariate
+ (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
+ female ies 1, then modmaxcovj=1.
+ */
+ } /* end for loop on individuals i */
+ printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
+ fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
+ cptcode=modmaxcovj;
+ /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
+ /*for (i=0; i<=cptcode; i++) {*/
+ for (j=modmincovj; j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */
+ printf("Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
+ fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
+ if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */
+ if( j != -1){
+ ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th
+ covariate for which somebody answered excluding
+ undefined. Usually 2: 0 and 1. */
+ }
+ ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th
+ covariate for which somebody answered including
+ undefined. Usually 3: -1, 0 and 1. */
+ }
+ /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for
+ * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
+ } /* Ndum[-1] number of undefined modalities */
+
+ /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
+ /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7.
+ If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125;
+ modmincovj=3; modmaxcovj = 7;
+ There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3;
+ which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10;
+ defining two dummy variables: variables V1_1 and V1_2.
+ nbcode[Tvar[j]][ij]=k;
+ nbcode[Tvar[j]][1]=0;
+ nbcode[Tvar[j]][2]=1;
+ nbcode[Tvar[j]][3]=2;
+ To be continued (not working yet).
+ */
+ ij=0; /* ij is similar to i but can jump over null modalities */
+ for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
+ if (Ndum[i] == 0) { /* If nobody responded to this modality k */
+ break;
+ }
+ ij++;
+ nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/
+ cptcode = ij; /* New max modality for covar j */
+ } /* end of loop on modality i=-1 to 1 or more */
+ break;
+ case 1: /* Testing on varying covariate, could be simple and
+ * should look at waves or product of fixed *
+ * varying. No time to test -1, assuming 0 and 1 only */
+ ij=0;
+ for(i=0; i<=1;i++){
+ nbcode[Tvar[k]][++ij]=i;
}
- ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
- covariate for which somebody answered including
- undefined. Usually 3: -1, 0 and 1. */
- }
- /* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for
- * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
- } /* Ndum[-1] number of undefined modalities */
-
- /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
- /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7.
- If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125;
- modmincovj=3; modmaxcovj = 7;
- There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3;
- which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10;
- defining two dummy variables: variables V1_1 and V1_2.
- nbcode[Tvar[j]][ij]=k;
- nbcode[Tvar[j]][1]=0;
- nbcode[Tvar[j]][2]=1;
- nbcode[Tvar[j]][3]=2;
- To be continued (not working yet).
- */
- ij=0; /* ij is similar to i but can jump over null modalities */
- for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
- if (Ndum[i] == 0) { /* If nobody responded to this modality k */
break;
- }
- ij++;
- nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
- cptcode = ij; /* New max modality for covar j */
- } /* end of loop on modality i=-1 to 1 or more */
+ default:
+ break;
+ } /* end switch */
+ } /* end dummy test */
/* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
/* /\*recode from 0 *\/ */
} /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/
for (k=-1; k< maxncov; k++) Ndum[k]=0;
-
+ /* Look at fixed dummy (single or product) covariates to check empty modalities */
for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */
/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
- ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
- Ndum[ij]++; /* Might be supersed V1 + V1*age */
+ ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */
+ Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */
+ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, {2, 1, 1, 1, 2, 1, 1, 0, 0} */
} /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */
ij=0;
- for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+ /* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */
+ for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
- if((Ndum[i]!=0) && (i<=ncovcol)){
+ /* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */
+ if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */
+ /* If product not in single variable we don't print results */
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
- Tvaraff[++ij]=i; /*For printing (unclear) */
- }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){
- Tvaraff[++ij]=-10; /* Dont'n know how to treat quantitative variables yet */
- }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){
- Tvaraff[++ij]=i; /*For printing (unclear) */
- }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){
- Tvaraff[++ij]=-20; /* Dont'n know how to treat quantitative variables yet */
- }
+ ++ij;
+ Tvaraff[ij]=Tvar[k]; /*For printing */
+ Tmodelind[ij]=k;
+ if(Fixed[k]!=0)
+ anyvaryingduminmodel=1;
+ /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */
+ /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */
+ /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */
+ /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */
+ /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */
+ /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */
+ }
} /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
/* ij--; */
/* cptcoveff=ij; /\*Number of total covariates*\/ */
*cptcov=ij; /*Number of total real effective covariates: effective
* because they can be excluded from the model and real
- * if in the model but excluded because missing values*/
+ * if in the model but excluded because missing values, but how to get k from ij?*/
+ for(j=ij+1; j<= cptcovt; j++){
+ Tvaraff[j]=0;
+ Tmodelind[j]=0;
+ }
+ /* To be sorted */
+ ;
}
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
if(k==cptcoveff){
- fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' with line ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
- 6+(cpt-1), cpt );
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
+ 4+(cpt-1), cpt ); /* 4 or 6 ?*/
}else{
fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
kl++;
fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
- /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
- /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
- /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
- vlv= nbcode[Tvaraff[k]][lv];
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
}
-
+
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
-set ter svg size 640, 480\n \
-unset log y\n \
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
k=3;
for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
- if(j==1)
- fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
- else
- fprintf(ficgp,", '' ");
- l=(nlstate+ndeath)*(cpt-1) +j;
- fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
- /* for (i=2; i<= nlstate+ndeath ; i ++) */
- /* fprintf(ficgp,"+$%d",k+l+i-1); */
- fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
+ if(j==1)
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+ else
+ fprintf(ficgp,", '' ");
+ l=(nlstate+ndeath)*(cpt-1) +j;
+ fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
+ /* for (i=2; i<= nlstate+ndeath ; i ++) */
+ /* fprintf(ficgp,"+$%d",k+l+i-1); */
+ fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
} /* nlstate */
fprintf(ficgp,", '' ");
fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
- l=(nlstate+ndeath)*(cpt-1) +j;
- if(j < nlstate)
- fprintf(ficgp,"$%d +",k+l);
- else
- fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
+ l=(nlstate+ndeath)*(cpt-1) +j;
+ if(j < nlstate)
+ fprintf(ficgp,"$%d +",k+l);
+ else
+ fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
}
fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
} /* end covariate */
-
+
/* 6eme */
/* CV preval stable (period) for each covariate */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-
+
fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
- /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
- /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
- /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
- vlv= nbcode[Tvaraff[k]][lv];
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}
fprintf(ficgp,"\n#\n");
if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
}
-
+
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n \
-unset log y\n \
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
k=3; /* Offset */
for (i=1; i<= nlstate ; i ++){
- if(i==1)
- fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
- else
- fprintf(ficgp,", '' ");
- l=(nlstate+ndeath)*(i-1)+1;
- fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
- for (j=2; j<= nlstate ; j ++)
- fprintf(ficgp,"+$%d",k+l+j-1);
- fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
+ if(i==1)
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+ else
+ fprintf(ficgp,", '' ");
+ l=(nlstate+ndeath)*(i-1)+1;
+ fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+ for (j=2; j<= nlstate ; j ++)
+ fprintf(ficgp,"+$%d",k+l+j-1);
+ fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
} /* nlstate */
fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
} /* end covariate */
-
-
+
+
/* 7eme */
if(backcast == 1){
/* CV back preval stable (period) for each covariate */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
- fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
- /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
- /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
- vlv= nbcode[Tvaraff[k]][lv];
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
- }
- fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
-
- fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
- fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n \
-unset log y\n \
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
+
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
- k=3; /* Offset */
- for (i=1; i<= nlstate ; i ++){
- if(i==1)
- fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
- else
- fprintf(ficgp,", '' ");
- /* l=(nlstate+ndeath)*(i-1)+1; */
- l=(nlstate+ndeath)*(cpt-1)+1;
- /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
- /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
- fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
- /* for (j=2; j<= nlstate ; j ++) */
- /* fprintf(ficgp,"+$%d",k+l+j-1); */
- /* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
- fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
- } /* nlstate */
- fprintf(ficgp,"\nset out\n");
+ k=3; /* Offset */
+ for (i=1; i<= nlstate ; i ++){
+ if(i==1)
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
+ else
+ fprintf(ficgp,", '' ");
+ /* l=(nlstate+ndeath)*(i-1)+1; */
+ l=(nlstate+ndeath)*(cpt-1)+1;
+ /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
+ /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
+ fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
+ /* for (j=2; j<= nlstate ; j ++) */
+ /* fprintf(ficgp,"+$%d",k+l+j-1); */
+ /* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
+ fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
} /* end covariate */
} /* End if backcast */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
- fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
- /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
- /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
- /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
- vlv= nbcode[Tvaraff[k]][lv];
- fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
- }
- fprintf(ficgp,"\n#\n");
- if(invalidvarcomb[k1]){
- fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- continue;
- }
-
- fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
- fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
- fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
-set ter svg size 640, 480\n \
-unset log y\n \
+ fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
+
+ fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
- for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
- /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
- /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
- /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
- /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
- if(i==1){
- fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
- }else{
- fprintf(ficgp,",\\\n '' ");
- }
- if(cptcoveff ==0){ /* No covariate */
- ioffset=2; /* Age is in 2 */
- /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
- /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
- /*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
- /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
- fprintf(ficgp," u %d:(", ioffset);
- if(i==nlstate+1)
- fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ", \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
- else
- fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
- }else{ /* more than 2 covariates */
- if(cptcoveff ==1){
- ioffset=4; /* Age is in 4 */
- }else{
- ioffset=6; /* Age is in 6 */
- /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
- /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
- }
- fprintf(ficgp," u %d:(",ioffset);
- kl=0;
- strcpy(gplotcondition,"(");
- for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
- /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
- /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
- /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
- vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */
- kl++;
- sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
- kl++;
- if(k <cptcoveff && cptcoveff>1)
- sprintf(gplotcondition+strlen(gplotcondition)," && ");
- }
- strcpy(gplotcondition+strlen(gplotcondition),")");
- /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
- /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
- /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
- /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
- if(i==nlstate+1){
- fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
- }else{
- fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
- }
- } /* end if covariate */
- } /* nlstate */
- fprintf(ficgp,"\nset out\n");
+ for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ if(i==1){
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
+ }else{
+ fprintf(ficgp,",\\\n '' ");
+ }
+ if(cptcoveff ==0){ /* No covariate */
+ ioffset=2; /* Age is in 2 */
+ /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ /*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ fprintf(ficgp," u %d:(", ioffset);
+ if(i==nlstate+1)
+ fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ else
+ fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ }else{ /* more than 2 covariates */
+ if(cptcoveff ==1){
+ ioffset=4; /* Age is in 4 */
+ }else{
+ ioffset=6; /* Age is in 6 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ }
+ fprintf(ficgp," u %d:(",ioffset);
+ kl=0;
+ strcpy(gplotcondition,"(");
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */
+ kl++;
+ sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
+ kl++;
+ if(k <cptcoveff && cptcoveff>1)
+ sprintf(gplotcondition+strlen(gplotcondition)," && ");
+ }
+ strcpy(gplotcondition+strlen(gplotcondition),")");
+ /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+ /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
+ /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
+ /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
+ if(i==nlstate+1){
+ fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ }else{
+ fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ }
+ } /* end if covariate */
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
} /* end covariate */
} /* End if prevfcast */
-
-
+
+
/* proba elementaires */
fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
for(i=1,jk=1; i <=nlstate; i++){
fprintf(ficgp,"# initial state %d\n",i);
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
- fprintf(ficgp,"# current state %d\n",k);
- for(j=1; j <=ncovmodel; j++){
- fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
- jk++;
- }
- fprintf(ficgp,"\n");
+ fprintf(ficgp,"# current state %d\n",k);
+ for(j=1; j <=ncovmodel; j++){
+ fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
+ jk++;
+ }
+ fprintf(ficgp,"\n");
}
}
}
fprintf(ficgp,"##############\n#\n");
-
+
/*goto avoid;*/
fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
}
}
else
- fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); /* Valgrind bug nbcode */
+ fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); /* Valgrind bug nbcode */
}
}else{
i=i-ncovmodel;
if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
fprintf(ficgp," (1.");
}
-
+
if(ng != 1){
fprintf(ficgp,")/(1");
-
+
for(k1=1; k1 <=nlstate; k1++){
if(nagesqr==0)
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
if(jprojmean==0) jprojmean=1;
if(mprojmean==0) jprojmean=1;
- i1=cptcoveff;
+ i1=pow(2,cptcoveff);
if (cptcovn < 1){i1=1;}
fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
fprintf(ficresf,"#****** Routine prevforecast **\n");
-
+
/* if (h==(int)(YEARM*yearp)){ */
- for(cptcov=1, k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
- k=k+1;
- fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- fprintf(ficresf," yearproj age");
- for(j=1; j<=nlstate+ndeath;j++){
- for(i=1; i<=nlstate;i++)
- fprintf(ficresf," p%d%d",i,j);
- fprintf(ficresf," wp.%d",j);
- }
- for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
- fprintf(ficresf,"\n");
- fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);
- for (agec=fage; agec>=(ageminpar-1); agec--){
- nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
- nhstepm = nhstepm/hstepm;
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- oldm=oldms;savm=savms;
- hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
-
- for (h=0; h<=nhstepm; h++){
- if (h*hstepm/YEARM*stepm ==yearp) {
- fprintf(ficresf,"\n");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
- }
- for(j=1; j<=nlstate+ndeath;j++) {
- ppij=0.;
- for(i=1; i<=nlstate;i++) {
- if (mobilav==1)
- ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
- else {
- ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
- }
- if (h*hstepm/YEARM*stepm== yearp) {
- fprintf(ficresf," %.3f", p3mat[i][j][h]);
- }
- } /* end i */
- if (h*hstepm/YEARM*stepm==yearp) {
- fprintf(ficresf," %.3f", ppij);
- }
- }/* end j */
- } /* end h */
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- } /* end agec */
- } /* end yearp */
- } /* end cptcod */
- } /* end cptcov */
+ for(k=1;k<=i1;k++){
+ if(invalidvarcomb[k]){
+ printf("\nCombination (%d) projection ignored because no cases \n",k);
+ continue;
+ }
+ fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficresf," yearproj age");
+ for(j=1; j<=nlstate+ndeath;j++){
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresf," p%d%d",i,j);
+ fprintf(ficresf," wp.%d",j);
+ }
+ for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
+ fprintf(ficresf,"\n");
+ fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);
+ for (agec=fage; agec>=(ageminpar-1); agec--){
+ nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
+ nhstepm = nhstepm/hstepm;
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+
+ for (h=0; h<=nhstepm; h++){
+ if (h*hstepm/YEARM*stepm ==yearp) {
+ fprintf(ficresf,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
+ }
+ for(j=1; j<=nlstate+ndeath;j++) {
+ ppij=0.;
+ for(i=1; i<=nlstate;i++) {
+ if (mobilav==1)
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
+ else {
+ ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k];
+ }
+ if (h*hstepm/YEARM*stepm== yearp) {
+ fprintf(ficresf," %.3f", p3mat[i][j][h]);
+ }
+ } /* end i */
+ if (h*hstepm/YEARM*stepm==yearp) {
+ fprintf(ficresf," %.3f", ppij);
+ }
+ }/* end j */
+ } /* end h */
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ } /* end agec */
+ } /* end yearp */
+ } /* end k */
fclose(ficresf);
printf("End of Computing forecasting \n");
/* } */
/************** Forecasting *****not tested NB*************/
-void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
+/* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */
- int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
- int *popage;
- double calagedatem, agelim, kk1, kk2;
- double *popeffectif,*popcount;
- double ***p3mat,***tabpop,***tabpopprev;
- /* double ***mobaverage; */
- char filerespop[FILENAMELENGTH];
+/* int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; */
+/* int *popage; */
+/* double calagedatem, agelim, kk1, kk2; */
+/* double *popeffectif,*popcount; */
+/* double ***p3mat,***tabpop,***tabpopprev; */
+/* /\* double ***mobaverage; *\/ */
+/* char filerespop[FILENAMELENGTH]; */
- tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- agelim=AGESUP;
- calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
+/* tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+/* tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+/* agelim=AGESUP; */
+/* calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; */
- prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+/* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
- strcpy(filerespop,"POP_");
- strcat(filerespop,fileresu);
- if((ficrespop=fopen(filerespop,"w"))==NULL) {
- printf("Problem with forecast resultfile: %s\n", filerespop);
- fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);
- }
- printf("Computing forecasting: result on file '%s' \n", filerespop);
- fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop);
+/* strcpy(filerespop,"POP_"); */
+/* strcat(filerespop,fileresu); */
+/* if((ficrespop=fopen(filerespop,"w"))==NULL) { */
+/* printf("Problem with forecast resultfile: %s\n", filerespop); */
+/* fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop); */
+/* } */
+/* printf("Computing forecasting: result on file '%s' \n", filerespop); */
+/* fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop); */
- if (cptcoveff==0) ncodemax[cptcoveff]=1;
+/* if (cptcoveff==0) ncodemax[cptcoveff]=1; */
- /* if (mobilav!=0) { */
- /* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
- /* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ */
- /* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
- /* printf(" Error in movingaverage mobilav=%d\n",mobilav); */
- /* } */
- /* } */
+/* /\* if (mobilav!=0) { *\/ */
+/* /\* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
+/* /\* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ *\/ */
+/* /\* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); *\/ */
+/* /\* printf(" Error in movingaverage mobilav=%d\n",mobilav); *\/ */
+/* /\* } *\/ */
+/* /\* } *\/ */
- stepsize=(int) (stepm+YEARM-1)/YEARM;
- if (stepm<=12) stepsize=1;
+/* stepsize=(int) (stepm+YEARM-1)/YEARM; */
+/* if (stepm<=12) stepsize=1; */
- agelim=AGESUP;
+/* agelim=AGESUP; */
- hstepm=1;
- hstepm=hstepm/stepm;
+/* hstepm=1; */
+/* hstepm=hstepm/stepm; */
- if (popforecast==1) {
- if((ficpop=fopen(popfile,"r"))==NULL) {
- printf("Problem with population file : %s\n",popfile);exit(0);
- fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0);
- }
- popage=ivector(0,AGESUP);
- popeffectif=vector(0,AGESUP);
- popcount=vector(0,AGESUP);
+/* if (popforecast==1) { */
+/* if((ficpop=fopen(popfile,"r"))==NULL) { */
+/* printf("Problem with population file : %s\n",popfile);exit(0); */
+/* fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0); */
+/* } */
+/* popage=ivector(0,AGESUP); */
+/* popeffectif=vector(0,AGESUP); */
+/* popcount=vector(0,AGESUP); */
- i=1;
- while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1;
+/* i=1; */
+/* while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1; */
- imx=i;
- for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
- }
+/* imx=i; */
+/* for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i]; */
+/* } */
- for(cptcov=1,k=0;cptcov<=i2;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
- k=k+1;
- fprintf(ficrespop,"\n#******");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- fprintf(ficrespop,"******\n");
- fprintf(ficrespop,"# Age");
- for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespop," P.%d",j);
- if (popforecast==1) fprintf(ficrespop," [Population]");
+/* for(cptcov=1,k=0;cptcov<=i2;cptcov++){ */
+/* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
+/* k=k+1; */
+/* fprintf(ficrespop,"\n#******"); */
+/* for(j=1;j<=cptcoveff;j++) { */
+/* fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
+/* } */
+/* fprintf(ficrespop,"******\n"); */
+/* fprintf(ficrespop,"# Age"); */
+/* for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespop," P.%d",j); */
+/* if (popforecast==1) fprintf(ficrespop," [Population]"); */
- for (cpt=0; cpt<=0;cpt++) {
- fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
+/* for (cpt=0; cpt<=0;cpt++) { */
+/* fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); */
- for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
- nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
- nhstepm = nhstepm/hstepm;
+/* for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ */
+/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); */
+/* nhstepm = nhstepm/hstepm; */
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- oldm=oldms;savm=savms;
- hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+/* oldm=oldms;savm=savms; */
+/* hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */
- for (h=0; h<=nhstepm; h++){
- if (h==(int) (calagedatem+YEARM*cpt)) {
- fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm);
- }
- for(j=1; j<=nlstate+ndeath;j++) {
- kk1=0.;kk2=0;
- for(i=1; i<=nlstate;i++) {
- if (mobilav==1)
- kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
- else {
- kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
- }
- }
- if (h==(int)(calagedatem+12*cpt)){
- tabpop[(int)(agedeb)][j][cptcod]=kk1;
- /*fprintf(ficrespop," %.3f", kk1);
- if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/
- }
- }
- for(i=1; i<=nlstate;i++){
- kk1=0.;
- for(j=1; j<=nlstate;j++){
- kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];
- }
- tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)];
- }
+/* for (h=0; h<=nhstepm; h++){ */
+/* if (h==(int) (calagedatem+YEARM*cpt)) { */
+/* fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); */
+/* } */
+/* for(j=1; j<=nlstate+ndeath;j++) { */
+/* kk1=0.;kk2=0; */
+/* for(i=1; i<=nlstate;i++) { */
+/* if (mobilav==1) */
+/* kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; */
+/* else { */
+/* kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; */
+/* } */
+/* } */
+/* if (h==(int)(calagedatem+12*cpt)){ */
+/* tabpop[(int)(agedeb)][j][cptcod]=kk1; */
+/* /\*fprintf(ficrespop," %.3f", kk1); */
+/* if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*\/ */
+/* } */
+/* } */
+/* for(i=1; i<=nlstate;i++){ */
+/* kk1=0.; */
+/* for(j=1; j<=nlstate;j++){ */
+/* kk1= kk1+tabpop[(int)(agedeb)][j][cptcod]; */
+/* } */
+/* tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)]; */
+/* } */
- if (h==(int)(calagedatem+12*cpt))
- for(j=1; j<=nlstate;j++)
- fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
- }
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- }
- }
+/* if (h==(int)(calagedatem+12*cpt)) */
+/* for(j=1; j<=nlstate;j++) */
+/* fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]); */
+/* } */
+/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+/* } */
+/* } */
- /******/
+/* /\******\/ */
- for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
- fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
- for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
- nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
- nhstepm = nhstepm/hstepm;
+/* for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) { */
+/* fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); */
+/* for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ */
+/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); */
+/* nhstepm = nhstepm/hstepm; */
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- oldm=oldms;savm=savms;
- hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
- for (h=0; h<=nhstepm; h++){
- if (h==(int) (calagedatem+YEARM*cpt)) {
- fprintf(ficresf,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm);
- }
- for(j=1; j<=nlstate+ndeath;j++) {
- kk1=0.;kk2=0;
- for(i=1; i<=nlstate;i++) {
- kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];
- }
- if (h==(int)(calagedatem+12*cpt)) fprintf(ficresf," %15.2f", kk1);
- }
- }
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- }
- }
- }
- }
+/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+/* oldm=oldms;savm=savms; */
+/* hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */
+/* for (h=0; h<=nhstepm; h++){ */
+/* if (h==(int) (calagedatem+YEARM*cpt)) { */
+/* fprintf(ficresf,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); */
+/* } */
+/* for(j=1; j<=nlstate+ndeath;j++) { */
+/* kk1=0.;kk2=0; */
+/* for(i=1; i<=nlstate;i++) { */
+/* kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod]; */
+/* } */
+/* if (h==(int)(calagedatem+12*cpt)) fprintf(ficresf," %15.2f", kk1); */
+/* } */
+/* } */
+/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+/* } */
+/* } */
+/* } */
+/* } */
- /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+/* /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
- if (popforecast==1) {
- free_ivector(popage,0,AGESUP);
- free_vector(popeffectif,0,AGESUP);
- free_vector(popcount,0,AGESUP);
- }
- free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- fclose(ficrespop);
-} /* End of popforecast */
+/* if (popforecast==1) { */
+/* free_ivector(popage,0,AGESUP); */
+/* free_vector(popeffectif,0,AGESUP); */
+/* free_vector(popcount,0,AGESUP); */
+/* } */
+/* free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+/* free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+/* fclose(ficrespop); */
+/* } /\* End of popforecast *\/ */
int fileappend(FILE *fichier, char *optionfich)
{
*/
{
int i, j, k, ks;
- int j1, k1, k2;
+ int j1, k1, k2, k3, k4;
char modelsav[80];
char stra[80], strb[80], strc[80], strd[80],stre[80];
char *strpt;
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */
/*
* Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
- for(k=cptcovt; k>=1;k--) /**< Number of covariates not including constant and age, neither age*age*/
- Tvar[k]=0;
+ for(k=cptcovt; k>=1;k--){ /**< Number of covariates not including constant and age, neither age*age*/
+ Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;
+ }
cptcovage=0;
for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
Typevar[k]=2; /* 2 for double fixed dummy covariates */
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */
+ Tposprod[k]=k1; /* Tpsprod[3]=1, Tposprod[2]=5 */
Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
/* scanf("%d",i);*/
cutl(strd,strc,strb,'V');
- ks++; /**< Number of simple covariates*/
+ ks++; /**< Number of simple covariates dummy or quantitative, fixe or varying */
cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */
Tvar[k]=atoi(strd);
Typevar[k]=0; /* 0 for simple covariates */
k = 1 2 3 4 5 6 7 8 9
Tvar[k]= 5 4 3 1+1+2+1+1=6 5 2 7 1 5
Typevar[k]= 0 0 0 2 1 0 2 1 1
- Fixed[Tvar[k]]1 1 1 1 2 0 1 2 3
- Dummy[Tvar[k]]1 0 0 0 2 1 1 2 3
+ Fixed[k] 1 1 1 1 3 0 0 or 2 2 3
+ Dummy[k] 1 0 0 0 3 1 1 2 3
+ Tmodelind[combination of covar]=k;
*/
/* Dispatching between quantitative and time varying covariates */
/* If Tvar[k] >ncovcol it is a product */
/* Tvar[k] is the value n of Vn with n varying for 1 to nvcol, or p Vp=Vn*Vm for product */
/* Computing effective variables, ie used by the model, that is from the cptcovt variables */
+ printf("Model=%s\n\
+Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\
+Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
+Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
+ fprintf(ficlog,"Model=%s\n\
+Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\
+Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
+Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
+
for(k=1, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
if (Tvar[k] <=ncovcol && (Typevar[k]==0 || Typevar[k]==2)){ /* Simple or product fixed dummy covariatee */
- Fixed[Tvar[k]]= 0;
- Dummy[Tvar[k]]= 0;
+ Fixed[k]= 0;
+ Dummy[k]= 0;
ncoveff++;
}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/
- Fixed[Tvar[k]]= 0;
- Dummy[Tvar[k]]= 1;
+ Fixed[k]= 0;
+ Dummy[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;
+ Fixed[k]= 1;
+ Dummy[k]= 0;
ntveff++; /* Only simple time varying dummy variable */
- }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
- if( Typevar[k]==0){
- Fixed[Tvar[k]]= 1;
- Dummy[Tvar[k]]= 1;
+ }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
nqtveff++;/* Only simple time varying quantitative variable */
+ }else if (Typevar[k] == 1) { /* product with age */
+ if (Tvar[k] <=ncovcol ){ /* Simple or product fixed dummy covariatee */
+ Fixed[k]= 2;
+ Dummy[k]= 2;
+ /* ncoveff++; */
+ }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
+ Fixed[k]= 2;
+ Dummy[k]= 3;
+ /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */
+ }else if( Tvar[k] <=ncovcol+nqv+ntv ){
+ Fixed[k]= 3;
+ Dummy[k]= 2;
+ /* ntveff++; /\* Only simple time varying dummy variable *\/ */
+ }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 3;
+ Dummy[k]= 3;
+ /* 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]);
- }
+ }else if (Typevar[k] == 2) { /* product without age */
+ k1=Tposprod[k];
+ if(Tvard[k1][1] <=ncovcol){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 0; /* or 2 ?*/
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }
+ }else if(Tvard[k1][1] <=ncovcol+nqv){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 0; /* or 2 ?*/
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 0; /* or 2 ?*/
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }
+ }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }
+ }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[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("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
}
- printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[Tvar[k]],Dummy[Tvar[k]]);
- fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[Tvar[k]],Dummy[Tvar[k]]);
+ printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]);
+ fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]);
+ }
+ /* Searching for doublons in the model */
+ for(k1=1; k1<= cptcovt;k1++){
+ for(k2=1; k2 <k1;k2++){
+ if((Typevar[k1]==Typevar[k2]) && (Fixed[Tvar[k1]]==Fixed[Tvar[k2]]) && (Dummy[Tvar[k1]]==Dummy[Tvar[k2]] )){
+ if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
+ if(Tvar[k1]==Tvar[k2]){
+ printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+ fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+ return(1);
+ }
+ }else if (Typevar[k1] ==2){
+ k3=Tposprod[k1];
+ k4=Tposprod[k2];
+ if( ((Tvard[k3][1]== Tvard[k4][1])&&(Tvard[k3][2]== Tvard[k4][2])) || ((Tvard[k3][1]== Tvard[k4][2])&&(Tvard[k3][2]== Tvard[k4][1])) ){
+ printf("Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+ fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+ return(1);
+ }
+ }
+ }
+ }
}
- 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 */
printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
}
- printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
- fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+ printf("\nComputing period (stable) prevalence: result on file '%s' \n", filerespl);
+ fprintf(ficlog,"\nComputing period (stable) prevalence: result on file '%s' \n", filerespl);
pstamp(ficrespl);
fprintf(ficrespl,"# Period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl);
fprintf(ficrespl,"#Age ");
agebase=ageminpar;
agelim=agemaxpar;
- i1=pow(2,ncoveff);
+ /* i1=pow(2,ncoveff); */
+ i1=pow(2,cptcoveff); /* Number of dummy covariates */
if (cptcovn < 1){i1=1;}
for(k=1; k<=i1;k++){
fprintf(ficrespl,"#******");
printf("#******");
fprintf(ficlog,"#******");
- for(j=1;j<=nqfveff;j++) {
- fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ for(j=1;j<=cptcoveff ;j++) {/* all covariates */
+ fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficrespl,"******\n");
printf("******\n");
fprintf(ficlog,"******\n");
- if(invalidvarcomb[k]){
- printf("\nCombination (%d) ignored because no cases \n",k);
- fprintf(ficrespl,"#Combination (%d) ignored because no cases \n",k);
- fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k);
+ if(invalidvarcomb[k]){
+ printf("\nCombination (%d) ignored because no case \n",k);
+ fprintf(ficrespl,"#Combination (%d) ignored because no case \n",k);
+ fprintf(ficlog,"\nCombination (%d) ignored because no case \n",k);
continue;
- }
+ }
fprintf(ficrespl,"#Age ");
- for(j=1;j<=nqfveff;j++) {
+ for(j=1;j<=cptcoveff;j++) {
fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i);
fprintf(ficrespl,"Total Years_to_converge\n");
-
+
for (age=agebase; age<=agelim; age++){
/* for (age=agebase; age<=agebase; age++){ */
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
fprintf(ficrespl,"%.0f ",age );
- for(j=1;j<=nqfveff;j++)
- fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
tot=0.;
for(i=1; i<=nlstate;i++){
- tot += prlim[i][i];
- fprintf(ficrespl," %.5f", prlim[i][i]);
+ tot += prlim[i][i];
+ fprintf(ficrespl," %.5f", prlim[i][i]);
}
fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
} /* Age */
agelim=agemaxpar;
- i1=pow(2,nqfveff);
+ i1=pow(2,cptcoveff);
if (cptcovn < 1){i1=1;}
-
- for(k=1; k<=i1;k++){
- /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
- /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
- //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- /* k=k+1; */
- /* to clean */
+
+ for(k=1; k<=i1;k++){
//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
fprintf(ficresplb,"#******");
printf("#******");
fprintf(ficlog,"#******");
- for(j=1;j<=nqfveff;j++) {
+ for(j=1;j<=cptcoveff ;j++) {/* all covariates */
fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficresplb,"******\n");
printf("******\n");
fprintf(ficlog,"******\n");
- if(invalidvarcomb[k]){
- printf("\nCombination (%d) ignored because no cases \n",k);
- fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k);
- fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k);
- continue;
- }
+ if(invalidvarcomb[k]){
+ printf("\nCombination (%d) ignored because no cases \n",k);
+ fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k);
+ fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k);
+ continue;
+ }
fprintf(ficresplb,"#Age ");
- for(j=1;j<=nqfveff;j++) {
+ for(j=1;j<=cptcoveff;j++) {
fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
for(i=1; i<=nlstate;i++) fprintf(ficresplb," %d-%d ",i,i);
if(mobilavproj > 0){
/* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
/* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
- bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
+ bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
}else if (mobilavproj == 0){
- printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
- fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
- exit(1);
+ printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+ fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+ exit(1);
}else{
- /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
- bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
+ /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+ bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
}
fprintf(ficresplb,"%.0f ",age );
- for(j=1;j<=nqfveff;j++)
- fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
tot=0.;
for(i=1; i<=nlstate;i++){
- tot += bprlim[i][i];
- fprintf(ficresplb," %.5f", bprlim[i][i]);
+ tot += bprlim[i][i];
+ fprintf(ficresplb," %.5f", bprlim[i][i]);
}
fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
} /* Age */
/* hstepm=1; aff par mois*/
pstamp(ficrespij);
fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
- i1= pow(2,nqfveff);
+ i1= pow(2,cptcoveff);
/* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
/* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
/* k=k+1; */
- for (k=1; k <= (int) pow(2,nqfveff); k++){
+ for (k=1; k <= (int) pow(2,cptcoveff); k++){
fprintf(ficrespij,"\n#****** ");
- for(j=1;j<=nqfveff;j++)
+ for(j=1;j<=cptcoveff;j++)
fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrespij,"******\n");
/* hstepm=1; aff par mois*/
pstamp(ficrespijb);
fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
- i1= pow(2,nqfveff);
+ i1= pow(2,cptcoveff);
/* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
/* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
/* k=k+1; */
- for (k=1; k <= (int) pow(2,nqfveff); k++){
+ for (k=1; k <= (int) pow(2,cptcoveff); k++){
fprintf(ficrespijb,"\n#****** ");
- for(j=1;j<=nqfveff;j++)
+ for(j=1;j<=cptcoveff;j++)
fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrespijb,"******\n");
if(invalidvarcomb[k]){
ncovcol + k1
If already ncovcol=4 and model=V2+V1+V1*V4+age*V3
Tvar[3=V1*V4]=4+1 etc */
- Tprod=ivector(1,NCOVMAX); /* Gives the position of a product */
+ Tprod=ivector(1,NCOVMAX); /* Gives the k position of the k1 product */
+ Tposprod=ivector(1,NCOVMAX); /* Gives the k1 product from the k position */
/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3
if V2+V1+V1*V4+age*V3+V3*V2 TProd[k1=2]=5 (V3*V2)
+ Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2
*/
Tvaraff=ivector(1,NCOVMAX); /* Unclear */
Tvard=imatrix(1,NCOVMAX,1,2); /* n=Tvard[k1][1] and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm
4 covariates (3 plus signs)
Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
*/
-
+ Tmodelind=ivector(1,NCOVMAX);/** five the k model position of an
+ * individual dummy, fixed or varying:
+ * Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4,
+ * 3, 1, 0, 0, 0, 0, 0, 0},
+ * model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/
/* Main decodemodel */
cptcoveff=0;
if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */
tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
- }
-
- ncovcombmax=pow(2,cptcoveff);
- invalidvarcomb=ivector(1, ncovcombmax);
- for(i=1;i<ncovcombmax;i++)
- invalidvarcomb[i]=0;
-
+ }
+
+ ncovcombmax=pow(2,cptcoveff);
+ invalidvarcomb=ivector(1, ncovcombmax);
+ for(i=1;i<ncovcombmax;i++)
+ invalidvarcomb[i]=0;
+
/* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in
V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
/* 1 to ncodemax[j] which is the maximum value of this jth covariate */
-
+
/* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
/*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
/* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
/* Calculates basic frequencies. Computes observed prevalence at single age
and for any valid combination of covariates
and prints on file fileres'p'. */
- freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
- firstpass, lastpass, stepm, weightopt, model);
+ freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
+ firstpass, lastpass, stepm, weightopt, model);
fprintf(fichtm,"\n");
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
/*#include "hpijx.h"*/
hPijx(p, bage, fage);
fclose(ficrespij);
-
+
/* ncovcombmax= pow(2,cptcoveff); */
/*-------------- Variance of one-step probabilities---*/
k=1;
varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
-
+
/* Prevalence for each covariates in probs[age][status][cov] */
probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
for(i=1;i<=AGESUP;i++)
prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
if (mobilav!=0 ||mobilavproj !=0 ) {
mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
- for(i=1;i<=AGESUP;i++)
- for(j=1;j<=nlstate;j++)
- for(k=1;k<=ncovcombmax;k++)
- mobaverages[i][j][k]=0.;
+ for(i=1;i<=AGESUP;i++)
+ for(j=1;j<=nlstate;j++)
+ for(k=1;k<=ncovcombmax;k++)
+ mobaverages[i][j][k]=0.;
mobaverage=mobaverages;
if (mobilav!=0) {
- if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
- fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
- printf(" Error in movingaverage mobilav=%d\n",mobilav);
- }
+ if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
+ fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+ printf(" Error in movingaverage mobilav=%d\n",mobilav);
+ }
}
/* /\* Prevalence for each covariates in probs[age][status][cov] *\/ */
/* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
else if (mobilavproj !=0) {
- if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
- fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
- printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
- }
+ if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
+ fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
+ printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
+ }
}
}/* end if moving average */
-
+
/*---------- Forecasting ------------------*/
/*if((stepm == 1) && (strcmp(model,".")==0)){*/
if(prevfcast==1){
printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
- for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ for (k=1; k <= (int) pow(2,cptcoveff); k++){ /* For any combination of dummy covariates, fixed and varying */
fprintf(ficreseij,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficreseij,"******\n");
printf("done evsij\n");fflush(stdout);
fprintf(ficlog,"done evsij\n");fflush(ficlog);
- /*---------- Health expectancies and variances ------------*/
+ /*---------- State-specific expectancies and variances ------------*/
strcpy(filerest,"T_");
strcpy(fileresstde,"STDE_");
strcat(fileresstde,fileresu);
if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
- printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
- fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
+ printf("Problem with State specific Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
+ fprintf(ficlog,"Problem with State specific Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
}
- printf(" Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
- fprintf(ficlog," Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+ printf(" Computing State-specific Expectancies and standard errors: result on file '%s' \n", fileresstde);
+ fprintf(ficlog," Computing State-specific Expectancies and standard errors: result on file '%s' \n", fileresstde);
strcpy(filerescve,"CVE_");
strcat(filerescve,fileresu);
if((ficrescveij=fopen(filerescve,"w"))==NULL) {
- printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
- fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
+ printf("Problem with Covar. State-specific Exp. resultfile: %s\n", filerescve); exit(0);
+ fprintf(ficlog,"Problem with Covar. State-specific Exp. resultfile: %s\n", filerescve); exit(0);
}
- printf(" Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
- fprintf(ficlog," Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+ printf(" Computing Covar. of State-specific Expectancies: result on file '%s' \n", filerescve);
+ fprintf(ficlog," Computing Covar. of State-specific Expectancies: result on file '%s' \n", filerescve);
strcpy(fileresv,"V_");
strcat(fileresv,fileresu);
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
}
- printf(" Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(stdout);
- fprintf(ficlog," Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(ficlog);
+ printf(" Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout);
+ fprintf(ficlog," Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog);
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ printf("\n#****** ");
fprintf(ficrest,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficlog,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++){
+ printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
fprintf(ficrest,"******\n");
+ fprintf(ficlog,"******\n");
+ printf("******\n");
fprintf(ficresstdeij,"\n#****** ");
fprintf(ficrescveij,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficresstdeij,"******\n");
fprintf(ficrescveij,"******\n");
fprintf(ficresvij,"\n#****** ");
for(j=1;j<=cptcoveff;j++)
- fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficresvij,"******\n");
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
- printf(" cvevsij %d, ",k);
- fprintf(ficlog, " cvevsij %d, ",k);
+ printf(" cvevsij combination#=%d, ",k);
+ fprintf(ficlog, " cvevsij combination#=%d, ",k);
cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
printf(" end cvevsij \n ");
fprintf(ficlog, " end cvevsij \n ");
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
- oldm=oldms;savm=savms; /* ZZ Segmentation fault */
- cptcod= 0; /* To be deleted */
- printf("varevsij %d \n",vpopbased);
- fprintf(ficlog, "varevsij %d \n",vpopbased);
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
- fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
- if(vpopbased==1)
- fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
- else
- fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
- fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
- for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
- fprintf(ficrest,"\n");
- /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
- epj=vector(1,nlstate+1);
- printf("Computing age specific period (stable) prevalences in each health state \n");
- fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
- for(age=bage; age <=fage ;age++){
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
- if (vpopbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][k];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][k];
- }
- }
+ oldm=oldms;savm=savms; /* ZZ Segmentation fault */
+ cptcod= 0; /* To be deleted */
+ printf("varevsij vpopbased=%d \n",vpopbased);
+ fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased);
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
+ if(vpopbased==1)
+ fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
+ else
+ fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
+ fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
+ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+ fprintf(ficrest,"\n");
+ /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
+ epj=vector(1,nlstate+1);
+ printf("Computing age specific period (stable) prevalences in each health state \n");
+ fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
+ for(age=bage; age <=fage ;age++){
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
+ if (vpopbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][k];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][k];
+ }
+ }
- fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
- /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
- /* printf(" age %4.0f ",age); */
- for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
- for(i=1, epj[j]=0.;i <=nlstate;i++) {
- epj[j] += prlim[i][i]*eij[i][j][(int)age];
- /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
- /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
- }
- epj[nlstate+1] +=epj[j];
- }
- /* printf(" age %4.0f \n",age); */
+ fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+ /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
+ /* printf(" age %4.0f ",age); */
+ for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
+ for(i=1, epj[j]=0.;i <=nlstate;i++) {
+ epj[j] += prlim[i][i]*eij[i][j][(int)age];
+ /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
+ }
+ epj[nlstate+1] +=epj[j];
+ }
+ /* printf(" age %4.0f \n",age); */
- for(i=1, vepp=0.;i <=nlstate;i++)
- for(j=1;j <=nlstate;j++)
- vepp += vareij[i][j][(int)age];
- fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
- for(j=1;j <=nlstate;j++){
- fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
- }
- fprintf(ficrest,"\n");
- }
+ for(i=1, vepp=0.;i <=nlstate;i++)
+ for(j=1;j <=nlstate;j++)
+ vepp += vareij[i][j][(int)age];
+ fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+ for(j=1;j <=nlstate;j++){
+ fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+ }
+ fprintf(ficrest,"\n");
+ }
} /* End vpopbased */
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
/*}*/
} /* End k */
- free_vector(weight,1,n);
- free_imatrix(Tvard,1,NCOVMAX,1,2);
- free_imatrix(s,1,maxwav+1,1,n);
- free_matrix(anint,1,maxwav,1,n);
- free_matrix(mint,1,maxwav,1,n);
- free_ivector(cod,1,n);
- free_ivector(tab,1,NCOVMAX);
- fclose(ficresstdeij);
- fclose(ficrescveij);
- fclose(ficresvij);
- fclose(ficrest);
- printf("done Health expectancies\n");fflush(stdout);
- fprintf(ficlog,"done Health expectancies\n");fflush(ficlog);
- fclose(ficpar);
-
- /*------- Variance of period (stable) prevalence------*/
+ printf("done State-specific expectancies\n");fflush(stdout);
+ fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
+
+ /*------- Variance of period (stable) prevalence------*/
+
strcpy(fileresvpl,"VPL_");
strcat(fileresvpl,fileresu);
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
}
printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
-
+
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-
+
for (k=1; k <= (int) pow(2,cptcoveff); k++){
- fprintf(ficresvpl,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresvpl,"******\n");
+ fprintf(ficresvpl,"\n#****** ");
+ printf("\n#****** ");
+ fprintf(ficlog,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficresvpl,"******\n");
+ printf("******\n");
+ fprintf(ficlog,"******\n");
- varpl=matrix(1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
- free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
+ varpl=matrix(1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+ varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
+ free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
/*}*/
}
-
+
fclose(ficresvpl);
printf("done variance-covariance of period prevalence\n");fflush(stdout);
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
-
+
+ free_vector(weight,1,n);
+ free_imatrix(Tvard,1,NCOVMAX,1,2);
+ free_imatrix(s,1,maxwav+1,1,n);
+ free_matrix(anint,1,maxwav,1,n);
+ free_matrix(mint,1,maxwav,1,n);
+ free_ivector(cod,1,n);
+ free_ivector(tab,1,NCOVMAX);
+ fclose(ficresstdeij);
+ fclose(ficrescveij);
+ fclose(ficresvij);
+ fclose(ficrest);
+ fclose(ficpar);
+
+
/*---------- End : free ----------------*/
if (mobilav!=0 ||mobilavproj !=0)
free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
} /* mle==-3 arrives here for freeing */
- /* endfree:*/
- free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
- free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
- free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
- free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
- free_ma3x(cotvar,1,maxwav,1,ntv,1,n);
- free_matrix(coqvar,1,maxwav,1,n);
- free_matrix(covar,0,NCOVMAX,1,n);
- free_matrix(matcov,1,npar,1,npar);
- free_matrix(hess,1,npar,1,npar);
- /*free_vector(delti,1,npar);*/
- free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
- free_matrix(agev,1,maxwav,1,imx);
- free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
-
- free_ivector(ncodemax,1,NCOVMAX);
- free_ivector(ncodemaxwundef,1,NCOVMAX);
- free_ivector(Dummy,-1,NCOVMAX);
- free_ivector(Fixed,-1,NCOVMAX);
- free_ivector(Typevar,-1,NCOVMAX);
- free_ivector(Tvar,1,NCOVMAX);
- free_ivector(Tprod,1,NCOVMAX);
- free_ivector(Tvaraff,1,NCOVMAX);
- free_ivector(invalidvarcomb,1,ncovcombmax);
- free_ivector(Tage,1,NCOVMAX);
-
- free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
- /* free_imatrix(codtab,1,100,1,10); */
+ /* endfree:*/
+ free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
+ free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
+ free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
+ free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
+ free_ma3x(cotvar,1,maxwav,1,ntv,1,n);
+ free_matrix(coqvar,1,maxwav,1,n);
+ free_matrix(covar,0,NCOVMAX,1,n);
+ free_matrix(matcov,1,npar,1,npar);
+ free_matrix(hess,1,npar,1,npar);
+ /*free_vector(delti,1,npar);*/
+ free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
+ free_matrix(agev,1,maxwav,1,imx);
+ free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
+
+ free_ivector(ncodemax,1,NCOVMAX);
+ free_ivector(ncodemaxwundef,1,NCOVMAX);
+ free_ivector(Dummy,-1,NCOVMAX);
+ free_ivector(Fixed,-1,NCOVMAX);
+ free_ivector(Typevar,-1,NCOVMAX);
+ free_ivector(Tvar,1,NCOVMAX);
+ free_ivector(Tposprod,1,NCOVMAX);
+ free_ivector(Tprod,1,NCOVMAX);
+ free_ivector(Tvaraff,1,NCOVMAX);
+ free_ivector(invalidvarcomb,1,ncovcombmax);
+ free_ivector(Tage,1,NCOVMAX);
+ free_ivector(Tmodelind,1,NCOVMAX);
+
+ free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
+ /* free_imatrix(codtab,1,100,1,10); */
fflush(fichtm);
fflush(ficgp);
-
+
if((nberr >0) || (nbwarn>0)){
printf("End of Imach with %d errors and/or %d warnings. Please look at the log file for details.\n",nberr,nbwarn);
fprintf(ficlog,"End of Imach with %d errors and/or warnings %d. Please look at the log file for details.\n",nberr,nbwarn);
printf("Local time at start %s\nLocal time at end %s",strstart, strtend);
fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend);
printf("Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout));
-
+
printf("Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time));
fprintf(ficlog,"Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout));
fprintf(ficlog,"Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time));
fclose(ficgp);
fclose(ficlog);
/*------ End -----------*/
-
-
- printf("Before Current directory %s!\n",pathcd);
+
+
+ printf("Before Current directory %s!\n",pathcd);
#ifdef WIN32
- if (_chdir(pathcd) != 0)
- printf("Can't move to directory %s!\n",path);
- if(_getcwd(pathcd,MAXLINE) > 0)
+ if (_chdir(pathcd) != 0)
+ printf("Can't move to directory %s!\n",path);
+ if(_getcwd(pathcd,MAXLINE) > 0)
#else
- if(chdir(pathcd) != 0)
- printf("Can't move to directory %s!\n", path);
- if (getcwd(pathcd, MAXLINE) > 0)
+ if(chdir(pathcd) != 0)
+ printf("Can't move to directory %s!\n", path);
+ if (getcwd(pathcd, MAXLINE) > 0)
#endif
printf("Current directory %s!\n",pathcd);
/*strcat(plotcmd,CHARSEPARATOR);*/
sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot);
printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout);
-
+
if((outcmd=system(plotcmd)) != 0){
printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd);
printf("\n Trying if gnuplot resides on the same directory that IMaCh\n");
else if (z[0] == 'g') system(plotcmd);
else if (z[0] == 'q') exit(0);
}
- end:
+end:
while (z[0] != 'q') {
printf("\nType q for exiting: "); fflush(stdout);
scanf("%s",z);
git://henry.ined.fr / .git / commitdiff