begin-prev-date,...
open gnuplot file
open html file
- period (stable) prevalence
- for age prevalim()
- h Pij x
- variance of p varprob
+ period (stable) prevalence | pl_nom 1-1 2-2 etc by covariate
+ for age prevalim() | #****** V1=0 V2=1 V3=1 V4=0 ******
+ | 65 1 0 2 1 3 1 4 0 0.96326 0.03674
+ freexexit2 possible for memory heap.
+
+ h Pij x | pij_nom ficrestpij
+ # Cov Agex agex+h hpijx with i,j= 1-1 1-2 1-3 2-1 2-2 2-3
+ 1 85 85 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000
+ 1 85 86 0.68299 0.22291 0.09410 0.71093 0.00000 0.28907
+
+ 1 65 99 0.00364 0.00322 0.99314 0.00350 0.00310 0.99340
+ 1 65 100 0.00214 0.00204 0.99581 0.00206 0.00196 0.99597
+ variance of p one-step probabilities varprob | prob_nom ficresprob #One-step probabilities and stand. devi in ()
+ Standard deviation of one-step probabilities | probcor_nom ficresprobcor #One-step probabilities and correlation matrix
+ Matrix of variance covariance of one-step probabilities | probcov_nom ficresprobcov #One-step probabilities and covariance matrix
+
forecasting if prevfcast==1 prevforecast call prevalence()
health expectancies
Variance-covariance of DFLE
#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
+#define codtabm(h,k) 1 & (h-1) >> (k-1) ;
#define MAXN 20000
#define YEARM 12. /**< Number of months per year */
#define AGESUP 130
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
int nvar=0, nforce=0; /* Number of variables, number of forces */
-int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov=0; /* Number of covariates, of covariates with '*age' */
+/* Number of covariates model=V2+V1+ V3*age+V2*V4 */
+int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */
+int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */
+int cptcovs=0; /**< cptcovs number of simple covariates V2+V1 =2 */
+int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */
+int cptcovprodnoage=0; /**< Number of covariate products without age */
+int cptcoveff=0; /* Total number of covariates to vary for printing results */
+int cptcov=0; /* Working variable */
int npar=NPARMAX;
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between
* wave mi and wave mi+1 is not an exact multiple of stepm. */
double jmean=1; /* Mean space between 2 waves */
+double **matprod2(); /* test */
double **oldm, **newm, **savm; /* Working pointers to matrices */
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
/*FILE *fic ; */ /* Used in readdata only */
double *weight;
int **s; /* Status */
double *agedc;
-double **covar; /**< covar[i,j], value of jth covariate for individual i,
+double **covar; /**< covar[j,i], value of jth covariate for individual i,
* covar=matrix(0,NCOVMAX,1,n);
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+int *Ndum; /** Freq of modality (tricode */
int **codtab; /**< codtab=imatrix(1,100,1,10); */
int **Tvard, *Tprod, cptcovprod, *Tvaraff;
double *lsurv, *lpop, *tpop;
return s;
}
+char *cutl(char *blocc, char *alocc, char *in, char occ)
+{
+ /* cuts string in into blocc and alocc where blocc ends before first occurence of char 'occ'
+ and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2')
+ gives blocc="abcdef2ghi" and alocc="j".
+ If occ is not found blocc is null and alocc is equal to in. Returns blocc
+ */
+ char *s, *t, *bl;
+ t=in;s=in;
+ while ((*in != occ) && (*in != '\0')){
+ *alocc++ = *in++;
+ }
+ if( *in == occ){
+ *(alocc)='\0';
+ s=++in;
+ }
+
+ if (s == t) {/* occ not found */
+ *(alocc-(in-s))='\0';
+ in=s;
+ }
+ while ( *in != '\0'){
+ *blocc++ = *in++;
+ }
+
+ *blocc='\0';
+ return t;
+}
char *cutv(char *blocc, char *alocc, char *in, char occ)
{
/* cuts string in into blocc and alocc where blocc ends before last occurence of char 'occ'
for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
return m;
- /* print *(*(m+1)+70) or print m[1][70]; print m+1 or print &(m[1])
+ /* print *(*(m+1)+70) or print m[1][70]; print m+1 or print &(m[1]) or &(m[1][0])
+m[i] = address of ith row of the table. &(m[i]) is its value which is another adress
+that of m[i][0]. In order to get the value p m[i][0] but it is unitialized.
*/
}
int i, ii,j,k;
double min, max, maxmin, maxmax,sumnew=0.;
- double **matprod2();
+ /* double **matprod2(); */ /* test */
double **out, cov[NCOVMAX+1], **pmij();
double **newm;
double agefin, delaymax=50 ; /* Max number of years to converge */
for (k=1; k<=cptcovn;k++) {
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
- /* printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/
+ /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
}
- for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
- for (k=1; k<=cptcovprod;k++)
- cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+ /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */
+ /* cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; */
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
+ /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+ /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */
savm=oldm;
sumnew=0;
for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
prlim[i][j]= newm[i][j]/(1-sumnew);
+ /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/
max=FMAX(max,prlim[i][j]);
min=FMIN(min,prlim[i][j]);
}
}
}
-
-/* for(ii=1; ii<= nlstate+ndeath; ii++){ */
-/* for(jj=1; jj<= nlstate+ndeath; jj++){ */
-/* printf("ddd %lf ",ps[ii][jj]); */
-/* } */
-/* printf("\n "); */
-/* } */
-/* printf("\n ");printf("%lf ",cov[2]); */
- /*
+
+ /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+ /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
+ /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+ /* } */
+ /* printf("\n "); */
+ /* } */
+ /* printf("\n ");printf("%lf ",cov[2]);*/
+ /*
for(i=1; i<= npar; i++) printf("%f ",x[i]);
goto end;*/
return ps;
/**************** Product of 2 matrices ******************/
-double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)
+double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b)
{
/* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */
/* in, b, out are matrice of pointers which should have been initialized
before: only the contents of out is modified. The function returns
a pointer to pointers identical to out */
- long i, j, k;
+ int i, j, k;
for(i=nrl; i<= nrh; i++)
- for(k=ncolol; k<=ncoloh; k++)
- for(j=ncl,out[i][k]=0.; j<=nch; j++)
- out[i][k] +=in[i][j]*b[j][k];
-
+ for(k=ncolol; k<=ncoloh; k++){
+ out[i][k]=0.;
+ for(j=ncl; j<=nch; j++)
+ out[i][k] +=in[i][j]*b[j][k];
+ }
return out;
}
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
for (k=1; k<=cptcovage;k++)
cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
- for (k=1; k<=cptcovprod;k++)
+ for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
to be observed in j being in i according to the model.
*/
- for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */
+ cov[2+k]=covar[Tvar[k]][i];
+ }
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]
is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]
has been calculated etc */
for (kk=1; kk<=cptcovage;kk++) {
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
}
+ /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */
+ /* 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */
savm=oldm;
oldm=newm;
} /* end mult */
fx=func(x);
for (i=1;i<=npar;i++) p2[i]=x[i];
- for(l=0 ; l <=lmax; l++){
+ for(l=0 ; l <=lmax; l++){ /* Enlarging the zone around the Maximum */
l1=pow(10,l);
delts=delt;
for(k=1 ; k <kmax; k=k+1){
delt = delta*(l1*k);
p2[theta]=x[theta] +delt;
- k1=func(p2)-fx;
+ k1=func(p2)-fx; /* Might be negative if too close to the theoretical maximum */
p2[theta]=x[theta]-delt;
k2=func(p2)-fx;
/*res= (k1-2.0*fx+k2)/delt/delt; */
first=1;
- for(k1=1; k1<=j ; k1++){ /* Loop on covariates */
- for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */
- j1++;
+ /* for(k1=1; k1<=j ; k1++){ /* Loop on covariates */
+ /* for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */
+ /* j1++;
+*/
+ for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/
for (i=-5; i<=nlstate+ndeath; i++)
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
for (z1=1; z1<=cptcoveff; z1++)
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){
+ /* Tests if the value of each of the covariates of i is equal to filter j1 */
bool=0;
- printf("bool=%d i=%d, z1=%d, i1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n",
- bool,i,z1, i1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1],
- j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);
+ /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n",
+ bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1],
+ j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/
}
}
/*}*/
}
}
- }
+ } /* end i */
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
pstamp(ficresp);
printf("Others in log...\n");
fprintf(ficlog,"\n");
}
- }
+ /*}*/
}
dateintmean=dateintsum/k2cpt;
double pos,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;
/* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
j1=0;
- j=cptcoveff;
+ /*j=cptcoveff;*/
if (cptcovn<1) {j=1;ncodemax[1]=1;}
- for(k1=1; k1<=j;k1++){
- for(i1=1; i1<=ncodemax[k1];i1++){
- j1++;
+ first=1;
+ for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
+ /*for(i1=1; i1<=ncodemax[k1];i1++){
+ j1++;*/
for (i=1; i<=nlstate; i++)
for(m=iagemin; m <= iagemax+3; m++)
}
}
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
- printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\n",jk,i,j1,probs[i][jk][j1]);
+ } 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 */
- } /* end i1 */
- } /* end k1 */
+ /*} *//* end i1 */
+ } /* end j1 */
/* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
/*free_vector(pp,1,nlstate);*/
}
/*********** Tricode ****************************/
-void tricode(int *Tvar, int **nbcode, int imx)
+void tricode(int *Tvar, int **nbcode, int imx, int *Ndum)
{
/**< Uses cptcovn+2*cptcovprod as the number of covariates */
/* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1
/* Boring subroutine which should only output nbcode[Tvar[j]][k]
- /* nbcode[Tvar[j][1]=
+ * Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2)
+ /* nbcode[Tvar[j]][1]=
*/
- int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
+ int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
int modmaxcovj=0; /* Modality max of covariates j */
+ int cptcode=0; /* Modality max of covariates j */
+ int modmincovj=0; /* Modality min of covariates j */
+
+
cptcoveff=0;
- for (k=0; k < maxncov; k++) Ndum[k]=0;
+ for (k=-1; k < maxncov; k++) Ndum[k]=0;
for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
- for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate j */
- for (i=1; i<=imx; i++) { /*reads the data file to get the maximum value of the
+ /* Loop on covariates without age and products */
+ for (j=1; j<=(cptcovs); j++) { /* model V1 + V2*age+ V3 + V3*V4 : V1 + V3 = 2 only */
+ for (i=1; i<=imx; i++) { /* Lopp on individuals: reads the data file to get the maximum value of the
modality of this covariate Vj*/
- ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Finds for covariate j, n=Tvar[j] of Vn . ij is the
+ 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);*/
- if (ij > modmaxcovj) modmaxcovj=ij;
/* 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 is 1, then modmaxcovj=1.*/
}
+ printf(" 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<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each modality of model-cov j */
- if( Ndum[i] != 0 )
- ncodemax[j]++;
- /* Number of modalities of the j th covariate. In fact
- ncodemax[j]=2 (dichotom. variables only) but it could be more for
- historical reasons */
+ /*for (i=0; i<=cptcode; i++) {*/
+ for (i=modmincovj; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */
+ printf("Frequencies of covariates %d V%d %d\n", j, Tvar[j], Ndum[i]);
+ if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */
+ ncodemax[j]++; /* ncodemax[j]= Number of non-null modalities of the j th covariate. */
+ }
+ /* 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 */
- ij=1;
- for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 for dichotomous */
- for (k=0; k<= modmaxcovj; k++) { /* k=-1 ? NCOVMAX*//* maxncov or modmaxcovj */
+ /* For covariate j, modalities could be 1, 2, 3, 4. If Ndum[2]=0 ncodemax[j] is not 4 but 3 */
+ /* If 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 (or 2 if 27 is too few) : ncodemax[j]=3;
+ which will be coded 0, 1, 2 which in binary on 3-1 digits are 0=00 1=01, 2=10; defining two dummy
+ 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;
+ */
+ ij=1; /* ij is similar to i but can jumps over null modalities */
+ for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 */
+ for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */
+ /*recode from 0 */
if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */
nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode.
k is a modality. If we have model=V1+V1*sex
} /* end of loop on modality */
} /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/
- for (k=0; k< maxncov; k++) Ndum[k]=0;
+ for (k=-1; k< maxncov; k++) Ndum[k]=0;
- for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and 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]++;
- }
+ for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and 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]++;
+ }
ij=1;
- for (i=1; i<= maxncov; 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) */
+ /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
if((Ndum[i]!=0) && (i<=ncovcol)){
- Tvaraff[ij]=i; /*For printing */
+ /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
+ Tvaraff[ij]=i; /*For printing (unclear) */
ij++;
- }
+ }else
+ Tvaraff[ij]=0;
}
ij--;
cptcoveff=ij; /*Number of total covariates*/
+
}
+
/*********** Health Expectancies ****************/
void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
free_vector(gmp,nlstate+1,nlstate+ndeath);
free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
- fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65");
+ fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240");
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
- fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 2 ",subdirf(fileresprobmorprev));
- fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l lt 3 ",subdirf(fileresprobmorprev));
- fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 3 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);
int i, j=0, i1, k1, l1, t, tj;
int k2, l2, j1, z1;
int k=0,l, cptcode;
- int first=1, first1;
+ int first=1, first1, first2;
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
double **dnewm,**doldm;
double *xp;
double *gp, *gm;
double **gradg, **trgradg;
double **mu;
- double age,agelim, cov[NCOVMAX];
+ double age,agelim, cov[NCOVMAX+1];
double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
int theta;
char fileresprob[FILENAMELENGTH];
char fileresprobcov[FILENAMELENGTH];
char fileresprobcor[FILENAMELENGTH];
-
double ***varpij;
strcpy(fileresprob,"prob");
To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.<br> \n");
cov[1]=1;
- tj=cptcoveff;
+ /* tj=cptcoveff; */
+ tj = (int) pow(2,cptcoveff);
if (cptcovn<1) {tj=1;ncodemax[1]=1;}
j1=0;
- for(t=1; t<=tj;t++){
- for(i1=1; i1<=ncodemax[t];i1++){
- j1++;
+ for(j1=1; j1<=tj;j1++){
+ /*for(i1=1; i1<=ncodemax[t];i1++){ */
+ /*j1++;*/
if (cptcovn>0) {
fprintf(ficresprob, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
fprintf(ficresprobcor, "**********\n#");
}
+ gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
+ trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
+ gp=vector(1,(nlstate)*(nlstate+ndeath));
+ gm=vector(1,(nlstate)*(nlstate+ndeath));
for (age=bage; age<=fage; age ++){
cov[2]=age;
for (k=1; k<=cptcovn;k++) {
- cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
+ cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4
+ * 1 1 1 1 1
+ * 2 2 1 1 1
+ * 3 1 2 1 1
+ */
+ /* nbcode[1][1]=0 nbcode[1][2]=1;*/
}
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
for (k=1; k<=cptcovprod;k++)
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
- gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
- trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
- gp=vector(1,(nlstate)*(nlstate+ndeath));
- gm=vector(1,(nlstate)*(nlstate+ndeath));
for(theta=1; theta <=npar; theta++){
for(i=1; i<=npar; i++)
matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
- free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
- free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
- free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
- free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
pmij(pmmij,cov,ncovmodel,x,nlstate);
i=0;
for (k=1; k<=(nlstate);k++){
for (l=1; l<=(nlstate+ndeath);l++){
- i=i++;
+ i++;
fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
for (j=1; j<=i;j++){
+ /* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */
fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
}
}
}/* end of loop for state */
} /* end of loop for age */
-
+ free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
+ free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
+ free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+ free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+
/* Confidence intervalle of pij */
/*
fprintf(ficgp,"\nunset parametric;unset label");
*/
/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
- first1=1;
+ first1=1;first2=2;
for (k2=1; k2<=(nlstate);k2++){
for (l2=1; l2<=(nlstate+ndeath);l2++){
if(l2==k2) continue;
lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
if ((lc2 <0) || (lc1 <0) ){
- printf("Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Continuing by making them positive: WRONG RESULTS.\n", lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);
- fprintf(ficlog,"Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e\n", lc1, lc2, v1, v2, cv12);fflush(ficlog);
- lc1=fabs(lc1);
- lc2=fabs(lc2);
+ if(first2==1){
+ first1=0;
+ printf("Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);
+ }
+ fprintf(ficlog,"Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog);
+ /* lc1=fabs(lc1); */ /* If we want to have them positive */
+ /* lc2=fabs(lc2); */
}
/* Eigen vectors */
first=0;
fprintf(ficgp,"\nset parametric;unset label");
fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
- fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+ fprintf(ficgp,"\nset ter png small size 320, 240");
fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
:<a href=\"%s%d%1d%1d-%1d%1d.png\">\
%s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,\
} /* k12 */
} /*l1 */
}/* k1 */
- } /* loop covariates */
+ /* } /* loop covariates */
}
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
- m=cptcoveff;
+ m=pow(2,cptcoveff);
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
jj1=0;
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}
/* Pij */
- fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d1.png\">%s%d1.png</a><br> \
-<img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d_1.png\">%s%d_1.png</a><br> \
+<img src=\"%s%d_1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
/* Quasi-incidences */
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
- before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d2.png\">%s%d2.png</a><br> \
-<img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d_2.png\">%s%d_2.png</a><br> \
+<img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
/* Period (stable) prevalence in each health state */
for(cpt=1; cpt<nlstate;cpt++){
- fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d%d.png\">%s%d%d.png</a><br> \
-<img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+ fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \
+<img src=\"%s%d_%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
}
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
fflush(fichtm);
fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
- m=cptcoveff;
+ m=pow(2,cptcoveff);
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
jj1=0;
}
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): %s%d%d.png <br>\
-<img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
+prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.png <br>\
+<img src=\"%s%d_%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
}
fprintf(fichtm,"\n<br>- Total life expectancy by age and \
health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
strcpy(optfileres,"vpl");
/* 1eme*/
for (cpt=1; cpt<= nlstate ; cpt ++) {
- for (k1=1; k1<= m ; k1 ++) {
- fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
- fprintf(ficgp,"\n#set out \"v%s%d%d.png\" \n",optionfilefiname,cpt,k1);
+ for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+ fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
+ fprintf(ficgp,"\n#set out \"v%s%d_%d.png\" \n",optionfilefiname,cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \n\
set ylabel \"Probability\" \n\
-set ter png small\n\
-set size 0.65,0.65\n\
+set ter png small size 320, 240\n\
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
}
- fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 1,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+ fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
}
- fprintf(ficgp,"\" t\"95\%% CI\" w l lt 2,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+ fprintf(ficgp,"\" t\"95\%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
}
- fprintf(ficgp,"\" t\"\" w l lt 2,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 3",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
+ fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
}
}
/*2 eme*/
for (k1=1; k1<= m ; k1 ++) {
fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),k1);
- fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
+ fprintf(ficgp,"set ylabel \"Years\" \nset ter png small size 320, 240\nplot [%.f:%.f] ",ageminpar,fage);
for (i=1; i<= nlstate+1 ; i ++) {
k=2*i;
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
}
- fprintf(ficgp,"\" t\"\" w l lt 1,");
+ fprintf(ficgp,"\" t\"\" w l lt 0,");
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
for (j=1; j<= nlstate+1 ; j ++) {
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
}
- if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 1");
- else fprintf(ficgp,"\" t\"\" w l lt 1,");
+ if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+ else fprintf(ficgp,"\" t\"\" w l lt 0,");
}
}
/* k=2+nlstate*(2*cpt-2); */
k=2+(nlstate+1)*(cpt-1);
fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
- fprintf(ficgp,"set ter png small\n\
-set size 0.65,0.65\n\
+ fprintf(ficgp,"set ter png small size 320, 240\n\
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),k1-1,k1-1,k,cpt);
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt<=nlstate ; cpt ++) {
k=3;
- fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
+ fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter png small\nset size 0.65,0.65\n\
+set ter png small size 320, 240\n\
unset log y\n\
plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,subdirf2(fileres,"pij"),k1,k+cpt+1,k+1);
}
}
}
-
+ /*goto avoid;*/
for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
for(jk=1; jk <=m; jk++) {
- fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng);
+ fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng);
if (ng==2)
fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
else
fprintf(ficgp,"\nset title \"Probability\"\n");
- fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar);
+ fprintf(ficgp,"\nset ter png small size 320, 240\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar);
i=1;
for(k2=1; k2<=nlstate; k2++) {
k3=i;
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
ij=1;/* To be checked else nbcode[0][0] wrong */
for(j=3; j <=ncovmodel; j++) {
- if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
- ij++;
- }
- else
+ /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /\* Bug valgrind *\/ */
+ /* /\*fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);*\/ */
+ /* ij++; */
+ /* } */
+ /* else */
fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
}
fprintf(ficgp,")/(1");
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
ij=1;
for(j=3; j <=ncovmodel; j++){
- if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
- fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
- ij++;
- }
- else
+ /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { */
+ /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); */
+ /* ij++; */
+ /* } */
+ /* else */
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
}
fprintf(ficgp,")");
} /* end k2 */
} /* end jk */
} /* end ng */
+ avoid:
fflush(ficgp);
} /* end gnuplot */
strcpy(optfileres,"vpl");
fprintf(ficgp,"set out \"graphmort.png\"\n ");
fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n ");
- fprintf(ficgp, "set ter png small\n set log y\n");
- fprintf(ficgp, "set size 0.65,0.65\n");
+ fprintf(ficgp, "set ter png small size 320, 240\n set log y\n");
+ /* fprintf(ficgp, "set size 0.65,0.65\n"); */
fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp);
}
cutv(stra, strb,line,' ');
if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
}
- else if(iout=sscanf(strb,"%s.") != 0){
+ else if(iout=sscanf(strb,"%s.",dummy) != 0){
month=99;
year=9999;
}else{
cutv(stra, strb,line,' ');
if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
}
- else if(iout=sscanf(strb,"%s.") != 0){
+ else if(iout=sscanf(strb,"%s.", dummy) != 0){
month=99;
year=9999;
}else{
+}
+void removespace(char *str) {
+ char *p1 = str, *p2 = str;
+ do
+ while (*p2 == ' ')
+ p2++;
+ while (*p1++ = *p2++);
}
-int decodemodel ( char model[], int lastobs)
+int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:
+ * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age
+ * - cptcovt total number of covariates of the model nbocc(+)+1 = 8
+ * - cptcovn or number of covariates k of the models excluding age*products =6
+ * - cptcovage number of covariates with age*products =2
+ * - cptcovs number of simple covariates
+ * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
+ * which is a new column after the 9 (ncovcol) variables.
+ * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual
+ * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
+ * Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
+ * - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
+ */
{
- int i, j, k;
+ int i, j, k, ks;
int i1, j1, k1, k2;
char modelsav[80];
- char stra[80], strb[80], strc[80], strd[80],stre[80];
+ char stra[80], strb[80], strc[80], strd[80],stre[80];
+ /*removespace(model);*/
if (strlen(model) >1){ /* If there is at least 1 covariate */
- j=0, j1=0, k1=1, k2=1;
- j=nbocc(model,'+'); /* j=Number of '+' */
- j1=nbocc(model,'*'); /* j1=Number of '*' */
- cptcovn=j+1; /* Number of covariates V1+V2*age+V3 =>(2 plus signs) + 1=3
- but the covariates which are product must be computed and stored. */
- cptcovprod=j1; /*Number of products V1*V2 +v3*age = 2 */
-
+ j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0;
+ j=nbocc(model,'+'); /**< j=Number of '+' */
+ j1=nbocc(model,'*'); /**< j1=Number of '*' */
+ cptcovs=j+1-j1; /**< Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2 */
+ cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/
+ /* including age products which are counted in cptcovage.
+ /* but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */
+ cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */
+ cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */
strcpy(modelsav,model);
if (strstr(model,"AGE") !=0){
printf("Error. AGE must be in lower case 'age' model=%s ",model);
return 1;
}
+ /* Design
+ * V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight
+ * < ncovcol=8 >
+ * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8
+ * k= 1 2 3 4 5 6 7 8
+ * cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
+ * covar[k,i], value of kth covariate if not including age for individual i:
+ * covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8)
+ * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[4]=3 Tvar[8]=8
+ * if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and
+ * Tage[++cptcovage]=k
+ * if products, new covar are created after ncovcol with k1
+ * Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11
+ * Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product
+ * Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8
+ * Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2];
+ * Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted
+ * V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11
+ * < ncovcol=8 >
+ * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2
+ * k= 1 2 3 4 5 6 7 8 9 10 11 12
+ * Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8
+ * p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6}
+ * p Tprod[1]@2={ 6, 5}
+ *p Tvard[1][1]@4= {7, 8, 5, 6}
+ * covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8
+ * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ *How to reorganize?
+ * Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age
+ * Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6}
+ * {2, 1, 4, 8, 5, 6, 3, 7}
+ * Struct []
+ */
+
/* This loop fills the array Tvar from the string 'model'.*/
/* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */
/* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */
/* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */
/* } */
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- for(k=cptcovn; k>=1;k--){
- cutv(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
- modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4
- */
+ /*
+ * 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 */
+ Tvar[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 '+'
+ modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
/*scanf("%d",i);*/
- if (strchr(strb,'*')) { /* Model includes a product V2+V1+V4+V3*age strb=V3*age */
- cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn: strb=V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
- if (strcmp(strc,"age")==0) { /* Vn*age */
+ if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
+ cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
+ if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
+ /* covar is not filled and then is empty */
cptcovprod--;
- cutv(strb,stre,strd,'V'); /* stre="V3" */
- Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=2 ; V1+V2*age Tvar[2]=2 */
+ cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
+ Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2 */
cptcovage++; /* Sums the number of covariates which include age as a product */
Tage[cptcovage]=k; /* Tage[1] = 4 */
/*printf("stre=%s ", stre);*/
} else if (strcmp(strd,"age")==0) { /* or age*Vn */
cptcovprod--;
- cutv(strb,stre,strc,'V');
+ cutl(stre,strb,strc,'V');
Tvar[k]=atoi(stre);
cptcovage++;
Tage[cptcovage]=k;
} else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/
/* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
- cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
- Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
+ cptcovn++;
+ cptcovprodnoage++;k1++;
+ cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+ Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
because this model-covariate is a construction we invent a new column
ncovcol + k1
If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
- cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */
+ 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 */
- Tvard[k1][1]=atoi(strc); /* m 1 for V1*/
- Tvard[k1][2]=atoi(stre); /* n 4 for V4*/
- Tvar[cptcovn+k2]=Tvard[k1][1]; /* Tvar[(cptcovn=4+k2=1)=5]= 1 (V1) */
- Tvar[cptcovn+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovn=4+(k2=1)+1)=6]= 4 (V4) */
+ Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
+ Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
+ k2=k2+2;
+ Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
+ Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
for (i=1; i<=lastobs;i++){
/* Computes the new covariate which is a product of
- covar[n][i]* covar[m][i] and stores it at ncovol+k1 */
+ covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
}
- k1++;
- k2=k2+2;
} /* End age is not in the model */
} /* End if model includes a product */
else { /* no more sum */
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
/* scanf("%d",i);*/
- cutv(strd,strc,strb,'V');
- Tvar[k]=atoi(strc);
+ cutl(strd,strc,strb,'V');
+ ks++; /**< Number of simple covariates */
+ cptcovn++;
+ Tvar[k]=atoi(strd);
}
strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
double kk1, kk2;
double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
double **ximort;
- char *alph[]={"a","a","b","c","d","e"}, str[4];
+ char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
int *dcwave;
char z[1]="c", occ;
ungetc(c,ficpar);
- covar=matrix(0,NCOVMAX,1,n);
+ covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
/* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
v1+v2*age+v2*v3 makes cptcovn = 3
*/
if (strlen(model)>1)
- cptcovn=nbocc(model,'+')+1;
- /* ncovprod */
- ncovmodel=2+cptcovn; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/
+ ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/
+ else
+ ncovmodel=2;
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */
npar= nforce*ncovmodel; /* Number of parameters like aij*/
matcov=matrix(1,npar,1,npar);
}
else{
- /* Read guess parameters */
+ /* Read guessed parameters */
/* Reads comments: lines beginning with '#' */
while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
}
fflush(ficlog);
+ /* Reads scales values */
p=param[1][1];
/* Reads comments: lines beginning with '#' */
}
fflush(ficlog);
+ /* Reads covariance matrix */
delti=delti3[1][1];
for(j=1; j <=npar; j++) matcov[i][j]=0.;
for(i=1; i <=npar; i++){
- fscanf(ficpar,"%s",&str);
+ fscanf(ficpar,"%s",str);
if(mle==1)
printf("%s",str);
fprintf(ficlog,"%s",str);
ncovcol + k1
If already ncovcol=4 and model=V2+V1+V1*V4+age*V3
Tvar[3=V1*V4]=4+1 etc */
- Tprod=ivector(1,15); /* Gives the position of a product */
+ Tprod=ivector(1,NCOVMAX); /* Gives the position of a product */
/* 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)
*/
- Tvaraff=ivector(1,15);
- Tvard=imatrix(1,15,1,2); /* n=Tvard[k1][1] and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm
+ 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
* For V3*V2 (in V2+V1+V1*V4+age*V3+V3*V2), V3*V2 position is 2nd.
* Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */
- Tage=ivector(1,15); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
+ Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
4 covariates (3 plus signs)
Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
*/
free_matrix(anint,1,maxwav,1,n);*/
free_vector(moisdc,1,n);
free_vector(andc,1,n);
-
-
+ /* */
+
wav=ivector(1,imx);
dh=imatrix(1,lastpass-firstpass+1,1,imx);
bh=imatrix(1,lastpass-firstpass+1,1,imx);
/* Concatenates waves */
concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm);
-
+ /* */
+
/* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
ncodemax[1]=1;
- if (cptcovn > 0) tricode(Tvar,nbcode,imx);
-
- codtab=imatrix(1,100,1,10); /**< codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1)
- */
+ Ndum =ivector(-1,NCOVMAX);
+ if (ncovmodel > 2)
+ tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
+
+ 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],codtab[100][10]);*/
h=0;
+
+
+ /*if (cptcovn > 0) */
+
+
m=pow(2,cptcoveff);
for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
* 16 2 2 2 1
*/
codtab[h][k]=j;
- codtab[h][Tvar[k]]=j;
+ /*codtab[h][Tvar[k]]=j;*/
printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);
}
}
printf("\n");
}
scanf("%d",i);*/
+
+ free_ivector(Ndum,-1,NCOVMAX);
+
+
/*------------ gnuplot -------------*/
strcpy(optionfilegnuplot,optionfilefiname);
ximort[i][j]=(i == j ? 1.0 : 0.0);
}
- p[1]=0.0268; p[NDIM]=0.083;
+ /*p[1]=0.0268; p[NDIM]=0.083;*/
/*printf("%lf %lf", p[1], p[2]);*/
- /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/
- /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
+ /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
+ /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
-
- strcpy(filerespl,"pl");
- strcat(filerespl,fileres);
- if((ficrespl=fopen(filerespl,"w"))==NULL) {
- printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
- fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
- }
- printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
- fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
- pstamp(ficrespl);
- fprintf(ficrespl,"# Period (stable) prevalence \n");
- fprintf(ficrespl,"#Age ");
- for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
- fprintf(ficrespl,"\n");
-
- prlim=matrix(1,nlstate,1,nlstate);
-
- agebase=ageminpar;
- agelim=agemaxpar;
- ftolpl=1.e-10;
- i1=pow(2,cptcoveff);
- if (cptcovn < 1){i1=1;}
-
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
- //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- /* to clean */
- //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]);
- fprintf(ficrespl,"\n#******");
- printf("\n#******");
- fprintf(ficlog,"\n#******");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- }
- fprintf(ficrespl,"******\n");
- printf("******\n");
- fprintf(ficlog,"******\n");
-
- for (age=agebase; age<=agelim; age++){
- prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
- fprintf(ficrespl,"%.0f ",age );
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- for(i=1; i<=nlstate;i++)
- fprintf(ficrespl," %.5f", prlim[i][i]);
- fprintf(ficrespl,"\n");
- } /* Age */
- /* was end of cptcod */
- } /* cptcov */
+#include "prevlim.h" /* Use ficrespl, ficlog */
fclose(ficrespl);
- /*------------- h Pij x at various ages ------------*/
-
- strcpy(filerespij,"pij"); strcat(filerespij,fileres);
- if((ficrespij=fopen(filerespij,"w"))==NULL) {
- printf("Problem with Pij resultfile: %s\n", filerespij);goto end;
- fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end;
- }
- printf("Computing pij: result on file '%s' \n", filerespij);
- fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);
-
- stepsize=(int) (stepm+YEARM-1)/YEARM;
- /*if (stepm<=24) stepsize=2;*/
-
- agelim=AGESUP;
- hstepm=stepsize*YEARM; /* Every year of age */
- hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
-
- /* 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 ");
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- fprintf(ficrespij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- fprintf(ficrespij,"******\n");
-
- for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
- nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
-
- /* nhstepm=nhstepm*YEARM; aff par mois*/
+#ifdef FREEEXIT2
+#include "freeexit2.h"
+#endif
- 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);
- fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespij," %1d-%1d",i,j);
- fprintf(ficrespij,"\n");
- for (h=0; h<=nhstepm; h++){
- fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespij," %.5f", p3mat[i][j][h]);
- fprintf(ficrespij,"\n");
- }
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- fprintf(ficrespij,"\n");
- }
- }
- }
+ /*------------- h Pij x at various ages ------------*/
+#include "hpijx.h"
+ fclose(ficrespij);
+ /*-------------- Variance of one-step probabilities---*/
+ k=1;
varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
- fclose(ficrespij);
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
for(i=1;i<=AGESUP;i++)
}
printf("Computing Health Expectancies: result on file '%s' \n", filerese);
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
+ /*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(ficreseij,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- }
+ /*}*/
}
fclose(ficreseij);
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- fprintf(ficrest,"\n#****** ");
+ /*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(ficrest,"\n#****** ");
for(j=1;j<=cptcoveff;j++)
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrest,"******\n");
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
+ /*
+ */
+ /* goto endfree; */
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
pstamp(ficrest);
+
+
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
- oldm=oldms;savm=savms;
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); 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 ");
+ oldm=oldms;savm=savms; /* Segmentation fault */
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart);
+ 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
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_vector(epj,1,nlstate+1);
- }
+ /*}*/
}
free_vector(weight,1,n);
- free_imatrix(Tvard,1,15,1,2);
+ 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);
}
printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- fprintf(ficresvpl,"\n#****** ");
+ /*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]][codtab[k][j]]);
fprintf(ficresvpl,"******\n");
oldm=oldms;savm=savms;
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k,strstart);
free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
- }
+ /*}*/
}
fclose(ficresvpl);
free_matrix(agev,1,maxwav,1,imx);
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
- free_ivector(ncodemax,1,8);
- free_ivector(Tvar,1,15);
- free_ivector(Tprod,1,15);
- free_ivector(Tvaraff,1,15);
- free_ivector(Tage,1,15);
+ free_ivector(ncodemax,1,NCOVMAX);
+ free_ivector(Tvar,1,NCOVMAX);
+ free_ivector(Tprod,1,NCOVMAX);
+ free_ivector(Tvaraff,1,NCOVMAX);
+ free_ivector(Tage,1,NCOVMAX);
free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
free_imatrix(codtab,1,100,1,10);
git://henry.ined.fr / .git / commitdiff