|
version 1.239, 2016/08/26 15:51:03
|
version 1.252, 2016/09/15 21:15:37
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.252 2016/09/15 21:15:37 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.251 2016/09/15 15:01:13 brouard |
| |
Summary: not working |
| |
|
| |
Revision 1.250 2016/09/08 16:07:27 brouard |
| |
Summary: continue |
| |
|
| |
Revision 1.249 2016/09/07 17:14:18 brouard |
| |
Summary: Starting values from frequencies |
| |
|
| |
Revision 1.248 2016/09/07 14:10:18 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.247 2016/09/02 11:11:21 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.246 2016/09/02 08:49:22 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.245 2016/09/02 07:25:01 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.244 2016/09/02 07:17:34 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.243 2016/09/02 06:45:35 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.242 2016/08/30 15:01:20 brouard |
| |
Summary: Fixing a lots |
| |
|
| |
Revision 1.241 2016/08/29 17:17:25 brouard |
| |
Summary: gnuplot problem in Back projection to fix |
| |
|
| |
Revision 1.240 2016/08/29 07:53:18 brouard |
| |
Summary: Better |
| |
|
| Revision 1.239 2016/08/26 15:51:03 brouard |
Revision 1.239 2016/08/26 15:51:03 brouard |
| Summary: Improvement in Powell output in order to copy and paste |
Summary: Improvement in Powell output in order to copy and paste |
| |
|
|
Line 92
|
Line 131
|
| Author: Nicolas Brouard |
Author: Nicolas Brouard |
| |
|
| Revision 1.210 2015/11/18 17:41:20 brouard |
Revision 1.210 2015/11/18 17:41:20 brouard |
| Summary: Start working on projected prevalences |
Summary: Start working on projected prevalences Revision 1.209 2015/11/17 22:12:03 brouard |
| |
|
| Revision 1.209 2015/11/17 22:12:03 brouard |
|
| Summary: Adding ftolpl parameter |
Summary: Adding ftolpl parameter |
| Author: N Brouard |
Author: N Brouard |
| |
|
|
Line 882 typedef struct {
|
Line 919 typedef struct {
|
| /* #include <libintl.h> */ |
/* #include <libintl.h> */ |
| /* #define _(String) gettext (String) */ |
/* #define _(String) gettext (String) */ |
| |
|
| #define MAXLINE 1024 /* Was 256. Overflow with 312 with 2 states and 4 covariates. Should be ok */ |
#define MAXLINE 2048 /* Was 256 and 1024. Overflow with 312 with 2 states and 4 covariates. Should be ok */ |
| |
|
| #define GNUPLOTPROGRAM "gnuplot" |
#define GNUPLOTPROGRAM "gnuplot" |
| /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ |
/*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ |
|
Line 1037 double dval;
|
Line 1074 double dval;
|
| #define FTOL 1.0e-10 |
#define FTOL 1.0e-10 |
| |
|
| #define NRANSI |
#define NRANSI |
| #define ITMAX 200 |
#define ITMAX 200 |
| |
#define ITPOWMAX 20 /* This is now multiplied by the number of parameters */ |
| |
|
| #define TOL 2.0e-4 |
#define TOL 2.0e-4 |
| |
|
|
Line 2066 void powell(double p[], double **xi, int
|
Line 2104 void powell(double p[], double **xi, int
|
| void linmin(double p[], double xi[], int n, double *fret, |
void linmin(double p[], double xi[], int n, double *fret, |
| double (*func)(double [])); |
double (*func)(double [])); |
| #else |
#else |
| void linmin(double p[], double xi[], int n, double *fret, |
void linmin(double p[], double xi[], int n, double *fret, |
| double (*func)(double []),int *flat); |
double (*func)(double []),int *flat); |
| #endif |
#endif |
| int i,ibig,j,jk,k; |
int i,ibig,j,jk,k; |
| double del,t,*pt,*ptt,*xit; |
double del,t,*pt,*ptt,*xit; |
|
Line 2107 void powell(double p[], double **xi, int
|
Line 2145 void powell(double p[], double **xi, int
|
| printf("\n#model= 1 + age "); |
printf("\n#model= 1 + age "); |
| fprintf(ficlog,"\n#model= 1 + age "); |
fprintf(ficlog,"\n#model= 1 + age "); |
| if(nagesqr==1){ |
if(nagesqr==1){ |
| printf(" + age*age ",Tvar[j]); |
printf(" + age*age "); |
| fprintf(ficlog," + age*age ",Tvar[j]); |
fprintf(ficlog," + age*age "); |
| } |
} |
| for(j=1;j <=ncovmodel-2;j++){ |
for(j=1;j <=ncovmodel-2;j++){ |
| if(Typevar[j]==0) { |
if(Typevar[j]==0) { |
|
Line 2131 void powell(double p[], double **xi, int
|
Line 2169 void powell(double p[], double **xi, int
|
| if (k != i) { |
if (k != i) { |
| printf("%d%d ",i,k); |
printf("%d%d ",i,k); |
| fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
| fprintf(ficres,"%1d%1d ",i,k); |
|
| for(j=1; j <=ncovmodel; j++){ |
for(j=1; j <=ncovmodel; j++){ |
| printf("%12.7f ",p[jk]); |
printf("%12.7f ",p[jk]); |
| fprintf(ficlog,"%12.7f ",p[jk]); |
fprintf(ficlog,"%12.7f ",p[jk]); |
| fprintf(ficres,"%12.7f ",p[jk]); |
|
| jk++; |
jk++; |
| } |
} |
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| fprintf(ficres,"\n"); |
|
| } |
} |
| } |
} |
| } |
} |
| if(*iter <=3){ |
if(*iter <=3 && *iter >1){ |
| tml = *localtime(&rcurr_time); |
tml = *localtime(&rcurr_time); |
| strcpy(strcurr,asctime(&tml)); |
strcpy(strcurr,asctime(&tml)); |
| rforecast_time=rcurr_time; |
rforecast_time=rcurr_time; |
| itmp = strlen(strcurr); |
itmp = strlen(strcurr); |
| if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ |
if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ |
| strcurr[itmp-1]='\0'; |
strcurr[itmp-1]='\0'; |
| printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); |
printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); |
| fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); |
fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); |
| for(niterf=10;niterf<=30;niterf+=10){ |
for(niterf=10;niterf<=30;niterf+=10){ |
| rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); |
rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); |
| forecast_time = *localtime(&rforecast_time); |
forecast_time = *localtime(&rforecast_time); |
| strcpy(strfor,asctime(&forecast_time)); |
strcpy(strfor,asctime(&forecast_time)); |
| itmp = strlen(strfor); |
itmp = strlen(strfor); |
| if(strfor[itmp-1]=='\n') |
if(strfor[itmp-1]=='\n') |
| strfor[itmp-1]='\0'; |
strfor[itmp-1]='\0'; |
| printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); |
printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); |
| fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); |
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); |
| } |
} |
| } |
} |
| for (i=1;i<=n;i++) { /* For each direction i */ |
for (i=1;i<=n;i++) { /* For each direction i */ |
|
Line 2215 void powell(double p[], double **xi, int
|
Line 2250 void powell(double p[], double **xi, int
|
| /* printf("\n"); */ |
/* printf("\n"); */ |
| /* fprintf(ficlog,"\n"); */ |
/* fprintf(ficlog,"\n"); */ |
| } |
} |
| if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */ |
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */ |
| |
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */ |
| /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ |
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ |
| /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */ |
/* By adding age*age in a model, the new -2LL should be lower and the difference follows a */ |
| /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */ |
/* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */ |
|
Line 2261 void powell(double p[], double **xi, int
|
Line 2297 void powell(double p[], double **xi, int
|
| free_vector(pt,1,n); |
free_vector(pt,1,n); |
| return; |
return; |
| } /* enough precision */ |
} /* enough precision */ |
| if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); |
if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations."); |
| for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ |
for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ |
| ptt[j]=2.0*p[j]-pt[j]; |
ptt[j]=2.0*p[j]-pt[j]; |
| xit[j]=p[j]-pt[j]; |
xit[j]=p[j]-pt[j]; |
|
Line 2541 Earliest age to start was %d-%d=%d, ncvl
|
Line 2577 Earliest age to start was %d-%d=%d, ncvl
|
| |
|
| /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */ |
/* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */ |
| /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */ |
/* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */ |
| double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij) |
double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij, int nres) |
| { |
{ |
| /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit |
/* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit |
| matrix by transitions matrix until convergence is reached with precision ftolpl */ |
matrix by transitions matrix until convergence is reached with precision ftolpl */ |
|
Line 2564 Earliest age to start was %d-%d=%d, ncvl
|
Line 2600 Earliest age to start was %d-%d=%d, ncvl
|
| /* If we start from prlim again, prlim tends to a constant matrix */ |
/* If we start from prlim again, prlim tends to a constant matrix */ |
| |
|
| int i, ii,j,k; |
int i, ii,j,k; |
| |
int first=0; |
| double *min, *max, *meandiff, maxmax,sumnew=0.; |
double *min, *max, *meandiff, maxmax,sumnew=0.; |
| /* double **matprod2(); */ /* test */ |
/* double **matprod2(); */ /* test */ |
| double **out, cov[NCOVMAX+1], **bmij(); |
double **out, cov[NCOVMAX+1], **bmij(); |
|
Line 2600 Earliest age to start was %d-%d=%d, ncvl
|
Line 2637 Earliest age to start was %d-%d=%d, ncvl
|
| cov[2]=agefin; |
cov[2]=agefin; |
| if(nagesqr==1) |
if(nagesqr==1) |
| cov[3]= agefin*agefin;; |
cov[3]= agefin*agefin;; |
| for (k=1; k<=cptcovn;k++) { |
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
| /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ |
/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ |
| cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
| /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */ |
/* printf("bprevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
| |
} |
| |
/* for (k=1; k<=cptcovn;k++) { */ |
| |
/* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ |
| |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ |
| |
/* /\* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); *\/ */ |
| |
/* } */ |
| |
for (k=1; k<=nsq;k++) { /* For single varying covariates only */ |
| |
/* Here comes the value of quantitative after renumbering k with single quantitative covariates */ |
| |
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; |
| |
/* printf("prevalim Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ |
| |
} |
| |
/* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; */ |
| |
/* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */ |
| |
/* /\* 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)]; */ |
| |
for (k=1; k<=cptcovage;k++){ /* For product with age */ |
| |
if(Dummy[Tvar[Tage[k]]]){ |
| |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
| |
} else{ |
| |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
| |
} |
| |
/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
| |
} |
| |
for (k=1; k<=cptcovprod;k++){ /* For product without age */ |
| |
/* printf("prevalim Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */ |
| |
if(Dummy[Tvard[k][1]==0]){ |
| |
if(Dummy[Tvard[k][2]==0]){ |
| |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; |
| |
}else{ |
| |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; |
| |
} |
| |
}else{ |
| |
if(Dummy[Tvard[k][2]==0]){ |
| |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; |
| |
}else{ |
| |
cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; |
| |
} |
| |
} |
| } |
} |
| for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; |
|
| for (k=1; k<=cptcovprod;k++) /* Useless */ |
|
| /* 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)]; |
|
| |
|
| /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ |
/*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 cov[4]=%lf \n",ij, cov[3],cov[4]);*/ |
|
Line 2655 Earliest age to start was %d-%d=%d, ncvl
|
Line 2726 Earliest age to start was %d-%d=%d, ncvl
|
| } |
} |
| } /* age loop */ |
} /* age loop */ |
| /* After some age loop it doesn't converge */ |
/* After some age loop it doesn't converge */ |
| printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ |
if(first){ |
| |
first=1; |
| |
printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\ |
| |
Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
| |
} |
| |
fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ |
| Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
| /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
| free_vector(min,1,nlstate); |
free_vector(min,1,nlstate); |
|
Line 2788 double **pmij(double **ps, double *cov,
|
Line 2864 double **pmij(double **ps, double *cov,
|
| /* }else */ |
/* }else */ |
| doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
| }else{ |
}else{ |
| printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); |
; |
| |
/* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */ |
| } |
} |
| } /*End ii */ |
} /*End ii */ |
| } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ |
} /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ |
|
Line 3171 double func( double *x)
|
Line 3248 double func( double *x)
|
| Then computes with function pmij which return a matrix p[i][j] giving the elementary probability |
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. |
to be observed in j being in i according to the model. |
| */ |
*/ |
| ioffset=2+nagesqr+cptcovage; |
ioffset=2+nagesqr ; |
| /* Fixed */ |
/* Fixed */ |
| for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */ |
for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */ |
| cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ |
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ |
|
Line 3190 double func( double *x)
|
Line 3267 double func( double *x)
|
| */ |
*/ |
| for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
| for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ |
for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ |
| cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; |
/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */ |
| |
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; |
| } |
} |
| for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
|
Line 3204 double func( double *x)
|
Line 3282 double func( double *x)
|
| if(nagesqr==1) |
if(nagesqr==1) |
| cov[3]= agexact*agexact; /* Should be changed here */ |
cov[3]= agexact*agexact; /* Should be changed here */ |
| for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
| |
if(!FixedV[Tvar[Tage[kk]]]) |
| cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ |
| |
else |
| |
cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; |
| } |
} |
| out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
| 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
|
Line 3486 double funcone( double *x)
|
Line 3567 double funcone( double *x)
|
| for(k=1; k<=nlstate; k++) ll[k]=0.; |
for(k=1; k<=nlstate; k++) ll[k]=0.; |
| ioffset=0; |
ioffset=0; |
| for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
| ioffset=2+nagesqr+cptcovage; |
/* ioffset=2+nagesqr+cptcovage; */ |
| |
ioffset=2+nagesqr; |
| /* Fixed */ |
/* Fixed */ |
| /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ |
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ |
| /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */ |
/* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */ |
|
Line 3513 double funcone( double *x)
|
Line 3595 double funcone( double *x)
|
| for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ |
for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ |
| /* Wave varying (but not age varying) */ |
/* Wave varying (but not age varying) */ |
| for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ |
for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ |
| cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; |
/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */ |
| } |
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; |
| |
} |
| /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */ |
/* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */ |
| /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ |
/* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ |
| /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ |
/* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ |
| /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */ |
/* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */ |
| /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */ |
/* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */ |
| /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ |
/* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ |
| /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */ |
/* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */ |
| /* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */ |
/* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */ |
| /* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */ |
/* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */ |
| /* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */ |
/* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */ |
| /* } */ |
/* } */ |
| for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
| savm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
| } |
} |
| |
|
| agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ |
agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ |
| ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ |
ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ |
| for(d=0; d<dh[mi][i]; d++){ /* Delay between two effective waves */ |
for(d=0; d<dh[mi][i]; d++){ /* Delay between two effective waves */ |
| /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] |
/* for(d=0; d<=0; d++){ /\* Delay between two effective waves Only one matrix to speed up*\/ */ |
| and mw[mi+1][i]. dh depends on stepm.*/ |
/*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] |
| newm=savm; |
and mw[mi+1][i]. dh depends on stepm.*/ |
| agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
newm=savm; |
| cov[2]=agexact; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; /* Here d is needed */ |
| if(nagesqr==1) |
cov[2]=agexact; |
| cov[3]= agexact*agexact; |
if(nagesqr==1) |
| for (kk=1; kk<=cptcovage;kk++) { |
cov[3]= agexact*agexact; |
| cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
| } |
if(!FixedV[Tvar[Tage[kk]]]) |
| /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */ |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
| /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
else |
| out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; |
| 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
} |
| /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */ |
/* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */ |
| /* 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */ |
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
| savm=oldm; |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
| oldm=newm; |
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 */ |
} /* end mult */ |
| |
|
| s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
|
Line 3566 double funcone( double *x)
|
Line 3653 double funcone( double *x)
|
| * is higher than the multiple of stepm and negative otherwise. |
* is higher than the multiple of stepm and negative otherwise. |
| */ |
*/ |
| if( s2 > nlstate && (mle <5) ){ /* Jackson */ |
if( s2 > nlstate && (mle <5) ){ /* Jackson */ |
| lli=log(out[s1][s2] - savm[s1][s2]); |
lli=log(out[s1][s2] - savm[s1][s2]); |
| } else if ( s2==-1 ) { /* alive */ |
} else if ( s2==-1 ) { /* alive */ |
| for (j=1,survp=0. ; j<=nlstate; j++) |
for (j=1,survp=0. ; j<=nlstate; j++) |
| survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
| lli= log(survp); |
lli= log(survp); |
| }else if (mle==1){ |
}else if (mle==1){ |
| lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
| } else if(mle==2){ |
} else if(mle==2){ |
| lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ |
lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ |
| } else if(mle==3){ /* exponential inter-extrapolation */ |
} else if(mle==3){ /* exponential inter-extrapolation */ |
| lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ |
lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ |
| } else if (mle==4){ /* mle=4 no inter-extrapolation */ |
} else if (mle==4){ /* mle=4 no inter-extrapolation */ |
| lli=log(out[s1][s2]); /* Original formula */ |
lli=log(out[s1][s2]); /* Original formula */ |
| } else{ /* mle=0 back to 1 */ |
} else{ /* mle=0 back to 1 */ |
| lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
| /*lli=log(out[s1][s2]); */ /* Original formula */ |
/*lli=log(out[s1][s2]); */ /* Original formula */ |
| } /* End of if */ |
} /* End of if */ |
| ipmx +=1; |
ipmx +=1; |
| sw += weight[i]; |
sw += weight[i]; |
| ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ |
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ |
| if(globpr){ |
if(globpr){ |
| fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ |
fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ |
| %11.6f %11.6f %11.6f ", \ |
%11.6f %11.6f %11.6f ", \ |
| num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, |
num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, |
| 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); |
2*weight[i]*lli,out[s1][s2],savm[s1][s2]); |
| for(k=1,llt=0.,l=0.; k<=nlstate; k++){ |
for(k=1,llt=0.,l=0.; k<=nlstate; k++){ |
| llt +=ll[k]*gipmx/gsw; |
llt +=ll[k]*gipmx/gsw; |
| fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); |
fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); |
| } |
} |
| fprintf(ficresilk," %10.6f\n", -llt); |
fprintf(ficresilk," %10.6f\n", -llt); |
| } |
} |
| } /* end of wave */ |
} /* end of wave */ |
| } /* end of individual */ |
} /* end of individual */ |
|
Line 3946 double hessij( double x[], double **hess
|
Line 4033 double hessij( double x[], double **hess
|
| kmax=kmax+10; |
kmax=kmax+10; |
| } |
} |
| if(kmax >=10 || firstime ==1){ |
if(kmax >=10 || firstime ==1){ |
| printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol); |
printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); |
| fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol); |
fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); |
| printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); |
printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); |
| fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); |
fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); |
| } |
} |
|
Line 4103 void pstamp(FILE *fichier)
|
Line 4190 void pstamp(FILE *fichier)
|
| } |
} |
| |
|
| /************ Frequencies ********************/ |
/************ Frequencies ********************/ |
| void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ |
void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ |
| int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \ |
int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \ |
| int firstpass, int lastpass, int stepm, int weightopt, char model[]) |
int firstpass, int lastpass, int stepm, int weightopt, char model[]) |
| { /* Some frequencies */ |
{ /* Some frequencies as well as proposing some starting values */ |
| |
|
| int i, m, jk, j1, bool, z1,j, k, iv; |
int i, m, jk, j1, bool, z1,j, k, iv, jj=0; |
| int iind=0, iage=0; |
int iind=0, iage=0; |
| int mi; /* Effective wave */ |
int mi; /* Effective wave */ |
| int first; |
int first; |
|
Line 4121 void freqsummary(char fileres[], int ia
|
Line 4208 void freqsummary(char fileres[], int ia
|
| double agebegin, ageend; |
double agebegin, ageend; |
| |
|
| pp=vector(1,nlstate); |
pp=vector(1,nlstate); |
| prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
| posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ |
posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ |
| pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ |
pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ |
| /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ |
/* prop=matrix(1,nlstate,iagemin,iagemax+3); */ |
|
Line 4135 void freqsummary(char fileres[], int ia
|
Line 4222 void freqsummary(char fileres[], int ia
|
| fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); |
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); |
| exit(0); |
exit(0); |
| } |
} |
| |
|
| strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); |
strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); |
| if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { |
if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { |
| printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
|
Line 4145 void freqsummary(char fileres[], int ia
|
Line 4232 void freqsummary(char fileres[], int ia
|
| } |
} |
| else{ |
else{ |
| fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
| <hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n \ |
| Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
| fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
| } |
} |
| fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm); |
fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm); |
| |
|
| strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); |
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); |
| if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { |
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { |
| printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
| fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
| fflush(ficlog); |
fflush(ficlog); |
| exit(70); |
exit(70); |
| } |
} else{ |
| else{ |
|
| fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
| <hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n \ |
| Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
| fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
| } |
} |
| fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr); |
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr); |
| |
|
| freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
| j1=0; |
j1=0; |
| |
|
| /* j=ncoveff; /\* Only fixed dummy covariates *\/ */ |
/* j=ncoveff; /\* Only fixed dummy covariates *\/ */ |
| j=cptcoveff; /* Only dummy covariates of the model */ |
j=cptcoveff; /* Only dummy covariates of the model */ |
| if (cptcovn<1) {j=1;ncodemax[1]=1;} |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
| |
|
| first=1; |
|
| |
|
| /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: |
/* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: |
| reference=low_education V1=0,V2=0 |
reference=low_education V1=0,V2=0 |
| med_educ V1=1 V2=0, |
med_educ V1=1 V2=0, |
| high_educ V1=0 V2=1 |
high_educ V1=0 V2=1 |
| Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff |
Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff |
| */ |
*/ |
| |
dateintsum=0; |
| |
k2cpt=0; |
| |
|
| 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 */ |
for (j = 0; j <= cptcoveff; j+=cptcoveff){ /* j= 0 constant model */ |
| posproptt=0.; |
first=1; |
| /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */ |
| scanf("%d", i);*/ |
posproptt=0.; |
| for (i=-5; i<=nlstate+ndeath; i++) |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
| for (jk=-5; jk<=nlstate+ndeath; jk++) |
scanf("%d", i);*/ |
| for(m=iagemin; m <= iagemax+3; m++) |
for (i=-5; i<=nlstate+ndeath; i++) |
| freq[i][jk][m]=0; |
for (jk=-5; jk<=nlstate+ndeath; jk++) |
| |
for(m=iagemin; m <= iagemax+3; m++) |
| for (i=1; i<=nlstate; i++) { |
freq[i][jk][m]=0; |
| for(m=iagemin; m <= iagemax+3; m++) |
|
| prop[i][m]=0; |
for (i=1; i<=nlstate; i++) { |
| posprop[i]=0; |
for(m=iagemin; m <= iagemax+3; m++) |
| pospropt[i]=0; |
prop[i][m]=0; |
| } |
posprop[i]=0; |
| /* for (z1=1; z1<= nqfveff; z1++) { */ |
pospropt[i]=0; |
| /* meanq[z1]+=0.; */ |
} |
| /* for(m=1;m<=lastpass;m++){ */ |
/* for (z1=1; z1<= nqfveff; z1++) { */ |
| /* meanqt[m][z1]=0.; */ |
/* meanq[z1]+=0.; */ |
| /* } */ |
/* for(m=1;m<=lastpass;m++){ */ |
| /* } */ |
/* meanqt[m][z1]=0.; */ |
| |
/* } */ |
| dateintsum=0; |
/* } */ |
| k2cpt=0; |
|
| /* For that combination of covariate j1, we count and print the frequencies in one pass */ |
/* dateintsum=0; */ |
| for (iind=1; iind<=imx; iind++) { /* For each individual iind */ |
/* k2cpt=0; */ |
| bool=1; |
|
| if(anyvaryingduminmodel==0){ /* If All fixed covariates */ |
/* For that combination of covariate j1, we count and print the frequencies in one pass */ |
| if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
for (iind=1; iind<=imx; iind++) { /* For each individual iind */ |
| /* for (z1=1; z1<= nqfveff; z1++) { */ |
bool=1; |
| /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ |
if(j !=0){ |
| /* } */ |
if(anyvaryingduminmodel==0){ /* If All fixed covariates */ |
| for (z1=1; z1<=cptcoveff; z1++) { |
if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
| /* if(Tvaraff[z1] ==-20){ */ |
/* for (z1=1; z1<= nqfveff; z1++) { */ |
| /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
/* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ |
| /* }else if(Tvaraff[z1] ==-10){ */ |
/* } */ |
| /* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ |
| /* }else */ |
/* if(Tvaraff[z1] ==-20){ */ |
| if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
/* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
| /* Tests if this individual iind responded to j1 (V4=1 V3=0) */ |
/* }else if(Tvaraff[z1] ==-10){ */ |
| bool=0; |
/* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
| /* 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", |
/* }else */ |
| bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), |
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */ |
| j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
/* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */ |
| /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ |
| } /* Onlyf fixed */ |
/* 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", |
| } /* end z1 */ |
bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), |
| } /* cptcovn > 0 */ |
j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
| } /* end any */ |
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
| if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */ |
} /* Onlyf fixed */ |
| /* for(m=firstpass; m<=lastpass; m++){ */ |
} /* end z1 */ |
| for(mi=1; mi<wav[iind];mi++){ /* For that wave */ |
} /* cptcovn > 0 */ |
| m=mw[mi][iind]; |
} /* end any */ |
| if(anyvaryingduminmodel==1){ /* Some are varying covariates */ |
}/* end j==0 */ |
| for (z1=1; z1<=cptcoveff; z1++) { |
if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */ |
| if( Fixed[Tmodelind[z1]]==1){ |
/* for(m=firstpass; m<=lastpass; m++){ */ |
| iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; |
for(mi=1; mi<wav[iind];mi++){ /* For that wave */ |
| if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */ |
m=mw[mi][iind]; |
| bool=0; |
if(j!=0){ |
| }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */ |
if(anyvaryingduminmodel==1){ /* Some are varying covariates */ |
| if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) { |
for (z1=1; z1<=cptcoveff; z1++) { |
| bool=0; |
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. If covariate's |
| |
value is -1, we don't select. It differs from the |
| |
constant and age model which counts them. */ |
| |
bool=0; /* not selected */ |
| |
}else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */ |
| |
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) { |
| |
bool=0; |
| |
} |
| |
} |
| } |
} |
| |
}/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ |
| |
} /* end j==0 */ |
| |
/* 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 */ |
| |
/* if((int)agev[m][iind] == 55) */ |
| |
/* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ |
| |
/* 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) && (j==0)) { |
| |
dateintsum=dateintsum+k2; /* on all covariates ?*/ |
| |
k2cpt++; |
| |
/* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ |
| } |
} |
| |
}else{ |
| |
bool=1; |
| |
}/* end bool 2 */ |
| |
} /* end m */ |
| |
} /* end bool */ |
| |
} /* end iind = 1 to imx */ |
| |
/* prop[s][age] is feeded for any initial and valid live state as well as |
| |
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ |
| |
|
| |
|
| |
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
| |
pstamp(ficresp); |
| |
if (cptcoveff>0 && j!=0){ |
| |
printf( "\n#********** Variable "); |
| |
fprintf(ficresp, "\n#********** Variable "); |
| |
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
| |
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); |
| |
fprintf(ficlog, "\n#********** Variable "); |
| |
for (z1=1; z1<=cptcoveff; z1++){ |
| |
if(!FixedV[Tvaraff[z1]]){ |
| |
printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresphtm, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresphtmfr, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficlog, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
}else{ |
| |
printf( "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
} |
| |
} |
| |
printf( "**********\n#"); |
| |
fprintf(ficresp, "**********\n#"); |
| |
fprintf(ficresphtm, "**********</h3>\n"); |
| |
fprintf(ficresphtmfr, "**********</h3>\n"); |
| |
fprintf(ficlog, "**********\n"); |
| |
} |
| |
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
| |
for(i=1; i<=nlstate;i++) { |
| |
fprintf(ficresp, " Age Prev(%d) N(%d) N ",i,i); |
| |
fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); |
| |
} |
| |
fprintf(ficresp, "\n"); |
| |
fprintf(ficresphtm, "\n"); |
| |
|
| |
/* Header of frequency table by age */ |
| |
fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">"); |
| |
fprintf(ficresphtmfr,"<th>Age</th> "); |
| |
for(jk=-1; jk <=nlstate+ndeath; jk++){ |
| |
for(m=-1; m <=nlstate+ndeath; m++){ |
| |
if(jk!=0 && m!=0) |
| |
fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m); |
| |
} |
| |
} |
| |
fprintf(ficresphtmfr, "\n"); |
| |
|
| |
/* For each age */ |
| |
for(iage=iagemin; iage <= iagemax+3; iage++){ |
| |
fprintf(ficresphtm,"<tr>"); |
| |
if(iage==iagemax+1){ |
| |
fprintf(ficlog,"1"); |
| |
fprintf(ficresphtmfr,"<tr><th>0</th> "); |
| |
}else if(iage==iagemax+2){ |
| |
fprintf(ficlog,"0"); |
| |
fprintf(ficresphtmfr,"<tr><th>Unknown</th> "); |
| |
}else if(iage==iagemax+3){ |
| |
fprintf(ficlog,"Total"); |
| |
fprintf(ficresphtmfr,"<tr><th>Total</th> "); |
| |
}else{ |
| |
if(first==1){ |
| |
first=0; |
| |
printf("See log file for details...\n"); |
| |
} |
| |
fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage); |
| |
fprintf(ficlog,"Age %d", iage); |
| |
} |
| |
for(jk=1; jk <=nlstate ; jk++){ |
| |
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) |
| |
pp[jk] += freq[jk][m][iage]; |
| |
} |
| |
for(jk=1; jk <=nlstate ; jk++){ |
| |
for(m=-1, pos=0; m <=0 ; m++) |
| |
pos += freq[jk][m][iage]; |
| |
if(pp[jk]>=1.e-10){ |
| |
if(first==1){ |
| |
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
| |
} |
| |
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
| |
}else{ |
| |
if(first==1) |
| |
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
| |
fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
| |
} |
| |
} |
| |
|
| |
for(jk=1; jk <=nlstate ; jk++){ |
| |
/* posprop[jk]=0; */ |
| |
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ |
| |
pp[jk] += freq[jk][m][iage]; |
| |
} /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ |
| |
|
| |
for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ |
| |
pos += pp[jk]; /* pos is the total number of transitions until this age */ |
| |
posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state |
| |
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
| |
pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state |
| |
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
| |
} |
| |
for(jk=1; jk <=nlstate ; jk++){ |
| |
if(pos>=1.e-5){ |
| |
if(first==1) |
| |
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
| |
fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
| |
}else{ |
| |
if(first==1) |
| |
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
| |
fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
| |
} |
| |
if( iage <= iagemax){ |
| |
if(pos>=1.e-5){ |
| |
fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
| |
fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
| |
/*probs[iage][jk][j1]= pp[jk]/pos;*/ |
| |
/*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ |
| |
} |
| |
else{ |
| |
fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); |
| |
fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta); |
| } |
} |
| }/* 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 */ |
pospropt[jk] +=posprop[jk]; |
| if(bool==1){ |
} /* end loop jk */ |
| /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] |
/* pospropt=0.; */ |
| and mw[mi+1][iind]. dh depends on stepm. */ |
for(jk=-1; jk <=nlstate+ndeath; jk++){ |
| agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ |
for(m=-1; m <=nlstate+ndeath; m++){ |
| ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */ |
if(freq[jk][m][iage] !=0 ) { /* minimizing output */ |
| if(m >=firstpass && m <=lastpass){ |
if(first==1){ |
| k2=anint[m][iind]+(mint[m][iind]/12.); |
printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); |
| /*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 */ |
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); |
| 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 */ |
if(jk!=0 && m!=0) |
| } /* end m */ |
fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]); |
| } /* end bool */ |
} |
| } /* end iind = 1 to imx */ |
} /* end loop jk */ |
| /* prop[s][age] is feeded for any initial and valid live state as well as |
posproptt=0.; |
| freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ |
for(jk=1; jk <=nlstate; jk++){ |
| |
posproptt += pospropt[jk]; |
| |
} |
| /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
fprintf(ficresphtmfr,"</tr>\n "); |
| pstamp(ficresp); |
if(iage <= iagemax){ |
| /* if (ncoveff>0) { */ |
fprintf(ficresp,"\n"); |
| if (cptcoveff>0) { |
fprintf(ficresphtm,"</tr>\n"); |
| fprintf(ficresp, "\n#********** Variable "); |
} |
| fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
if(first==1) |
| fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); |
printf("Others in log...\n"); |
| for (z1=1; z1<=cptcoveff; z1++){ |
fprintf(ficlog,"\n"); |
| fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
} /* end loop age iage */ |
| fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficresphtm,"<tr><th>Tot</th>"); |
| fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
|
| } |
|
| fprintf(ficresp, "**********\n#"); |
|
| fprintf(ficresphtm, "**********</h3>\n"); |
|
| fprintf(ficresphtmfr, "**********</h3>\n"); |
|
| fprintf(ficlog, "\n#********** Variable "); |
|
| for (z1=1; z1<=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\">"); |
|
| for(i=1; i<=nlstate;i++) { |
|
| fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); |
|
| fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); |
|
| } |
|
| fprintf(ficresp, "\n"); |
|
| fprintf(ficresphtm, "\n"); |
|
| |
|
| /* Header of frequency table by age */ |
|
| fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">"); |
|
| fprintf(ficresphtmfr,"<th>Age</th> "); |
|
| for(jk=-1; jk <=nlstate+ndeath; jk++){ |
|
| for(m=-1; m <=nlstate+ndeath; m++){ |
|
| if(jk!=0 && m!=0) |
|
| fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m); |
|
| } |
|
| } |
|
| fprintf(ficresphtmfr, "\n"); |
|
| |
|
| /* For each age */ |
|
| for(iage=iagemin; iage <= iagemax+3; iage++){ |
|
| fprintf(ficresphtm,"<tr>"); |
|
| if(iage==iagemax+1){ |
|
| fprintf(ficlog,"1"); |
|
| fprintf(ficresphtmfr,"<tr><th>0</th> "); |
|
| }else if(iage==iagemax+2){ |
|
| fprintf(ficlog,"0"); |
|
| fprintf(ficresphtmfr,"<tr><th>Unknown</th> "); |
|
| }else if(iage==iagemax+3){ |
|
| fprintf(ficlog,"Total"); |
|
| fprintf(ficresphtmfr,"<tr><th>Total</th> "); |
|
| }else{ |
|
| if(first==1){ |
|
| first=0; |
|
| printf("See log file for details...\n"); |
|
| } |
|
| fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage); |
|
| fprintf(ficlog,"Age %d", iage); |
|
| } |
|
| for(jk=1; jk <=nlstate ; jk++){ |
for(jk=1; jk <=nlstate ; jk++){ |
| for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) |
if(posproptt < 1.e-5){ |
| pp[jk] += freq[jk][m][iage]; |
fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt); |
| |
}else{ |
| |
fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt); |
| |
} |
| } |
} |
| for(jk=1; jk <=nlstate ; jk++){ |
fprintf(ficresphtm,"</tr>\n"); |
| for(m=-1, pos=0; m <=0 ; m++) |
fprintf(ficresphtm,"</table>\n"); |
| pos += freq[jk][m][iage]; |
fprintf(ficresphtmfr,"</table>\n"); |
| if(pp[jk]>=1.e-10){ |
if(posproptt < 1.e-5){ |
| if(first==1){ |
fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
| printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
| } |
fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); |
| fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
invalidvarcomb[j1]=1; |
| }else{ |
}else{ |
| if(first==1) |
fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1); |
| printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
invalidvarcomb[j1]=0; |
| fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
|
| } |
|
| } |
} |
| |
fprintf(ficresphtmfr,"</table>\n"); |
| for(jk=1; jk <=nlstate ; jk++){ |
fprintf(ficlog,"\n"); |
| /* posprop[jk]=0; */ |
if(j!=0){ |
| for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ |
printf("#Freqsummary: Starting values for combination j1=%d:\n", j1); |
| pp[jk] += freq[jk][m][iage]; |
for(i=1,jk=1; i <=nlstate; i++){ |
| } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ |
for(k=1; k <=(nlstate+ndeath); k++){ |
| |
if (k != i) { |
| for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ |
for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */ |
| pos += pp[jk]; /* pos is the total number of transitions until this age */ |
if(jj==1){ /* Constant case */ |
| posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state |
if(j1==1){ /* All dummy covariates to zero */ |
| from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */ |
| pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state |
freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */ |
| from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
printf("%d%d ",i,k); |
| |
fprintf(ficlog,"%d%d ",i,k); |
| |
printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3])); |
| |
fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
| |
pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); |
| |
} |
| |
}else if(jj==2 || nagesqr==1){ /* age or age*age parameter */ |
| |
; |
| |
}else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ |
| |
printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); |
| |
printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); |
| |
pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])); |
| |
printf("%d%d ",i,k); |
| |
fprintf(ficlog,"%d%d ",i,k); |
| |
printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4])); |
| |
}else{ /* Other cases, like quantitative fixed or varying covariates */ |
| |
; |
| |
} |
| |
/* printf("%12.7f )", param[i][jj][k]); */ |
| |
/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ |
| |
jk++; |
| |
} /* end jj */ |
| |
} /* end k!= i */ |
| |
} /* end k */ |
| |
} /* end i, jk */ |
| |
} /* end j !=0 */ |
| |
} /* end selected combination of covariate j1 */ |
| |
if(j==0){ /* We can estimate starting values from the occurences in each case */ |
| |
printf("#Freqsummary: Starting values for the constants:\n"); |
| |
fprintf(ficlog,"\n"); |
| |
for(i=1,jk=1; i <=nlstate; i++){ |
| |
for(k=1; k <=(nlstate+ndeath); k++){ |
| |
if (k != i) { |
| |
printf("%d%d ",i,k); |
| |
fprintf(ficlog,"%d%d ",i,k); |
| |
for(jj=1; jj <=ncovmodel; jj++){ |
| |
if(jj==1){ |
| |
printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
| |
fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
| |
} |
| |
/* printf("%12.7f )", param[i][jj][k]); */ |
| |
/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ |
| |
jk++; |
| |
} |
| |
printf("\n"); |
| |
fprintf(ficlog,"\n"); |
| |
} |
| |
} |
| } |
} |
| for(jk=1; jk <=nlstate ; jk++){ |
printf("#Freqsummary\n"); |
| if(pos>=1.e-5){ |
fprintf(ficlog,"\n"); |
| if(first==1) |
|
| printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
|
| fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
|
| }else{ |
|
| if(first==1) |
|
| printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
|
| fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
|
| } |
|
| if( iage <= iagemax){ |
|
| if(pos>=1.e-5){ |
|
| fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
|
| fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
|
| /*probs[iage][jk][j1]= pp[jk]/pos;*/ |
|
| /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ |
|
| } |
|
| else{ |
|
| fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); |
|
| fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta); |
|
| } |
|
| } |
|
| pospropt[jk] +=posprop[jk]; |
|
| } /* end loop jk */ |
|
| /* pospropt=0.; */ |
|
| for(jk=-1; jk <=nlstate+ndeath; jk++){ |
for(jk=-1; jk <=nlstate+ndeath; jk++){ |
| for(m=-1; m <=nlstate+ndeath; m++){ |
for(m=-1; m <=nlstate+ndeath; m++){ |
| if(freq[jk][m][iage] !=0 ) { /* minimizing output */ |
/* param[i]|j][k]= freq[jk][m][iagemax+3] */ |
| if(first==1){ |
printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); |
| printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); |
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); |
| } |
/* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */ |
| fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); |
/* printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ |
| } |
/* fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ |
| if(jk!=0 && m!=0) |
/* } */ |
| fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]); |
} |
| } |
|
| } /* end loop jk */ |
} /* end loop jk */ |
| posproptt=0.; |
|
| for(jk=1; jk <=nlstate; jk++){ |
printf("\n"); |
| posproptt += pospropt[jk]; |
|
| } |
|
| fprintf(ficresphtmfr,"</tr>\n "); |
|
| if(iage <= iagemax){ |
|
| fprintf(ficresp,"\n"); |
|
| fprintf(ficresphtm,"</tr>\n"); |
|
| } |
|
| if(first==1) |
|
| printf("Others in log...\n"); |
|
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| } /* end loop age iage */ |
} /* end j=0 */ |
| fprintf(ficresphtm,"<tr><th>Tot</th>"); |
} /* end j */ |
| for(jk=1; jk <=nlstate ; jk++){ |
|
| if(posproptt < 1.e-5){ |
if(mle == -2){ |
| fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt); |
for(i=1, jk=1; i <=nlstate; i++){ |
| }else{ |
for(j=1; j <=nlstate+ndeath; j++){ |
| fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt); |
if(j!=i){ |
| |
/*ca[0]= k+'a'-1;ca[1]='\0';*/ |
| |
printf("%1d%1d",i,j); |
| |
fprintf(ficparo,"%1d%1d",i,j); |
| |
for(k=1; k<=ncovmodel;k++){ |
| |
/* printf(" %lf",param[i][j][k]); */ |
| |
/* fprintf(ficparo," %lf",param[i][j][k]); */ |
| |
p[jk]=pstart[jk]; |
| |
printf(" %f ",pstart[jk]); |
| |
fprintf(ficparo," %f ",pstart[jk]); |
| |
jk++; |
| |
} |
| |
printf("\n"); |
| |
fprintf(ficparo,"\n"); |
| |
} |
| } |
} |
| } |
} |
| fprintf(ficresphtm,"</tr>\n"); |
} /* end mle=-2 */ |
| fprintf(ficresphtm,"</table>\n"); |
|
| fprintf(ficresphtmfr,"</table>\n"); |
|
| if(posproptt < 1.e-5){ |
|
| fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
|
| fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
|
| fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); |
|
| invalidvarcomb[j1]=1; |
|
| }else{ |
|
| fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1); |
|
| invalidvarcomb[j1]=0; |
|
| } |
|
| fprintf(ficresphtmfr,"</table>\n"); |
|
| } /* end selected combination of covariate j1 */ |
|
| dateintmean=dateintsum/k2cpt; |
dateintmean=dateintsum/k2cpt; |
| |
|
| fclose(ficresp); |
fclose(ficresp); |
| fclose(ficresphtm); |
fclose(ficresphtm); |
| fclose(ficresphtmfr); |
fclose(ficresphtmfr); |
| free_vector(meanq,1,nqfveff); |
free_vector(meanq,1,nqfveff); |
| free_matrix(meanqt,1,lastpass,1,nqtveff); |
free_matrix(meanqt,1,lastpass,1,nqtveff); |
| free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE); |
free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE); |
| free_vector(pospropt,1,nlstate); |
free_vector(pospropt,1,nlstate); |
| free_vector(posprop,1,nlstate); |
free_vector(posprop,1,nlstate); |
| free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE); |
free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+4+AGEMARGE); |
| free_vector(pp,1,nlstate); |
free_vector(pp,1,nlstate); |
| /* End of freqsummary */ |
/* End of freqsummary */ |
| } |
} |
|
Line 4483 void prevalence(double ***probs, double
|
Line 4692 void prevalence(double ***probs, double
|
| iagemin= (int) agemin; |
iagemin= (int) agemin; |
| iagemax= (int) agemax; |
iagemax= (int) agemax; |
| /*pp=vector(1,nlstate);*/ |
/*pp=vector(1,nlstate);*/ |
| prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
| /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ |
/* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ |
| j1=0; |
j1=0; |
| |
|
|
Line 4493 void prevalence(double ***probs, double
|
Line 4702 void prevalence(double ***probs, double
|
| first=1; |
first=1; |
| for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ |
for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ |
| for (i=1; i<=nlstate; i++) |
for (i=1; i<=nlstate; i++) |
| for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) |
for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) |
| prop[i][iage]=0.0; |
prop[i][iage]=0.0; |
| printf("Prevalence combination of varying and fixed dummies %d\n",j1); |
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," V%d=%d ",Tvaraff[j1],nbcode[Tvaraff[j1]][codtabm(k,j1)]); */ |
|
Line 4524 void prevalence(double ***probs, double
|
Line 4733 void prevalence(double ***probs, double
|
| if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ |
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]==0) agev[m][i]=iagemax+1; |
| if(agev[m][i]==1) agev[m][i]=iagemax+2; |
if(agev[m][i]==1) agev[m][i]=iagemax+2; |
| if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){ |
if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+4+AGEMARGE){ |
| printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); |
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); |
exit(1); |
| } |
} |
|
Line 4561 void prevalence(double ***probs, double
|
Line 4770 void prevalence(double ***probs, double
|
| |
|
| /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ |
/* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ |
| /*free_vector(pp,1,nlstate);*/ |
/*free_vector(pp,1,nlstate);*/ |
| free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
| } /* End of prevalence */ |
} /* End of prevalence */ |
| |
|
| /************* Waves Concatenation ***************/ |
/************* Waves Concatenation ***************/ |
|
Line 4782 void concatwav(int wav[], int **dh, int
|
Line 4991 void concatwav(int wav[], int **dh, int
|
| |
|
| /*********** Tricode ****************************/ |
/*********** Tricode ****************************/ |
| void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum) |
void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum) |
| { |
{ |
| /**< Uses cptcovn+2*cptcovprod as the number of covariates */ |
/**< 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 |
/* 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] |
* Boring subroutine which should only output nbcode[Tvar[j]][k] |
| * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable |
* Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable |
| * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually); |
* nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually); |
| */ |
*/ |
| |
|
| int 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 modmaxcovj=0; /* Modality max of covariates j */ |
| int cptcode=0; /* Modality max of covariates j */ |
int cptcode=0; /* Modality max of covariates j */ |
| int modmincovj=0; /* Modality min of covariates j */ |
int modmincovj=0; /* Modality min of covariates j */ |
| |
|
| |
|
| /* cptcoveff=0; */ |
/* cptcoveff=0; */ |
| /* *cptcov=0; */ |
/* *cptcov=0; */ |
| |
|
| for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ |
for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ |
| |
|
| /* Loop on covariates without age and products and no quantitative variable */ |
/* 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<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */ |
| for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ |
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; |
for (j=-1; (j < maxncov); j++) Ndum[j]=0; |
| if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
| switch(Fixed[k]) { |
switch(Fixed[k]) { |
| case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ |
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*/ |
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=(int)(covar[Tvar[k]][i]); |
| /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i |
/* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i |
| * If product of Vn*Vm, still boolean *: |
* If product of Vn*Vm, still boolean *: |
| * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables |
* 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 */ |
* 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 |
/* Finds for covariate j, n=Tvar[j] of Vn . ij is the |
| modality of the nth covariate of individual i. */ |
modality of the nth covariate of individual i. */ |
| if (ij > modmaxcovj) |
if (ij > modmaxcovj) |
| modmaxcovj=ij; |
modmaxcovj=ij; |
| else if (ij < modmincovj) |
else if (ij < modmincovj) |
| modmincovj=ij; |
modmincovj=ij; |
| if ((ij < -1) && (ij > NCOVMAX)){ |
if ((ij < -1) && (ij > NCOVMAX)){ |
| printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
| exit(1); |
exit(1); |
| }else |
}else |
| Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/ |
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 */ |
/* 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);*/ |
/*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 |
/* 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 |
(should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and |
| female ies 1, then modmaxcovj=1. |
female ies 1, then modmaxcovj=1. |
| */ |
*/ |
| } /* end for loop on individuals i */ |
} /* 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); |
printf(" Minimal and maximal values of %d th (fixed) 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); |
fprintf(ficlog," Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj); |
| cptcode=modmaxcovj; |
cptcode=modmaxcovj; |
| /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ |
/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ |
| /*for (i=0; i<=cptcode; i++) {*/ |
/*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 */ |
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]); |
printf("Frequencies of (fixed) covariate %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]); |
fprintf(ficlog, "Frequencies of (fixed) covariate %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( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */ |
| if( j != -1){ |
if( j != -1){ |
| ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th |
ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th |
| covariate for which somebody answered excluding |
covariate for which somebody answered excluding |
| undefined. Usually 2: 0 and 1. */ |
undefined. Usually 2: 0 and 1. */ |
| } |
} |
| ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th |
ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th |
| covariate for which somebody answered including |
covariate for which somebody answered including |
| undefined. Usually 3: -1, 0 and 1. */ |
undefined. Usually 3: -1, 0 and 1. */ |
| } /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for |
} /* 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 */ |
* historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ |
| } /* Ndum[-1] number of undefined modalities */ |
} /* Ndum[-1] number of undefined modalities */ |
| |
|
| /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ |
/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ |
| /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */ |
/* 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; */ |
/* 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; */ |
/* modmincovj=3; modmaxcovj = 7; */ |
| /* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */ |
/* 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; */ |
/* 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.*/ |
/* defining two dummy variables: variables V1_1 and V1_2.*/ |
| /* nbcode[Tvar[j]][ij]=k; */ |
/* nbcode[Tvar[j]][ij]=k; */ |
| /* nbcode[Tvar[j]][1]=0; */ |
/* nbcode[Tvar[j]][1]=0; */ |
| /* nbcode[Tvar[j]][2]=1; */ |
/* nbcode[Tvar[j]][2]=1; */ |
| /* nbcode[Tvar[j]][3]=2; */ |
/* nbcode[Tvar[j]][3]=2; */ |
| /* To be continued (not working yet). */ |
/* To be continued (not working yet). */ |
| ij=0; /* ij is similar to i but can jump over null modalities */ |
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*/ |
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 */ |
if (Ndum[i] == 0) { /* If nobody responded to this modality k */ |
| break; |
break; |
| } |
} |
| ij++; |
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*/ |
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 */ |
cptcode = ij; /* New max modality for covar j */ |
| } /* end of loop on modality i=-1 to 1 or more */ |
} /* end of loop on modality i=-1 to 1 or more */ |
| break; |
break; |
| case 1: /* Testing on varying covariate, could be simple and |
case 1: /* Testing on varying covariate, could be simple and |
| * should look at waves or product of fixed * |
* should look at waves or product of fixed * |
| * varying. No time to test -1, assuming 0 and 1 only */ |
* varying. No time to test -1, assuming 0 and 1 only */ |
| ij=0; |
ij=0; |
| for(i=0; i<=1;i++){ |
for(i=0; i<=1;i++){ |
| nbcode[Tvar[k]][++ij]=i; |
nbcode[Tvar[k]][++ij]=i; |
| } |
} |
| break; |
break; |
| default: |
default: |
| break; |
break; |
| } /* end switch */ |
} /* end switch */ |
| } /* end dummy test */ |
} /* end dummy test */ |
| |
|
| /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ |
/* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ |
| /* /\*recode from 0 *\/ */ |
/* /\*recode from 0 *\/ */ |
| /* k is a modality. If we have model=V1+V1*sex */ |
/* k is a modality. If we have model=V1+V1*sex */ |
| /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ |
/* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ |
| /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ |
/* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ |
| /* } */ |
/* } */ |
| /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ |
/* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ |
| /* if (ij > ncodemax[j]) { */ |
/* if (ij > ncodemax[j]) { */ |
| /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
/* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
| /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
/* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
| /* break; */ |
/* break; */ |
| /* } */ |
/* } */ |
| /* } /\* end of loop on modality k *\/ */ |
/* } /\* end of loop on modality k *\/ */ |
| } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ |
} /* 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; |
for (k=-1; k< maxncov; k++) Ndum[k]=0; |
| /* Look at fixed dummy (single or product) covariates to check empty modalities */ |
/* 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 */ |
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.*/ |
/* 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 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*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 */ |
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} */ |
/* 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} */ |
} /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */ |
| |
|
| ij=0; |
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) */ |
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]);*/ |
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ |
| /* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */ |
/* 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(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 */ |
/* If product not in single variable we don't print results */ |
| /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ |
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ |
| ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */ |
++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */ |
| Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/ |
Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/ |
| Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */ |
Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */ |
| TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ |
TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ |
| if(Fixed[k]!=0) |
if(Fixed[k]!=0) |
| anyvaryingduminmodel=1; |
anyvaryingduminmodel=1; |
| /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ |
/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ |
| /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ |
/* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ |
| /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ |
/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ |
| /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ |
/* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ |
| /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ |
/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ |
| /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ |
/* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ |
| } |
} |
| } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ |
} /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ |
| /* ij--; */ |
/* ij--; */ |
| /* cptcoveff=ij; /\*Number of total covariates*\/ */ |
/* cptcoveff=ij; /\*Number of total covariates*\/ */ |
| *cptcov=ij; /*Number of total real effective covariates: effective |
*cptcov=ij; /*Number of total real effective covariates: effective |
| * because they can be excluded from the model and real |
* because they can be excluded from the model and real |
| * if in the model but excluded because missing values, but how to get k from ij?*/ |
* if in the model but excluded because missing values, but how to get k from ij?*/ |
| for(j=ij+1; j<= cptcovt; j++){ |
for(j=ij+1; j<= cptcovt; j++){ |
| Tvaraff[j]=0; |
Tvaraff[j]=0; |
| Tmodelind[j]=0; |
Tmodelind[j]=0; |
| } |
} |
| for(j=ntveff+1; j<= cptcovt; j++){ |
for(j=ntveff+1; j<= cptcovt; j++){ |
| TmodelInvind[j]=0; |
TmodelInvind[j]=0; |
| } |
} |
| /* To be sorted */ |
/* To be sorted */ |
| ; |
; |
| } |
} |
| |
|
| |
|
| /*********** Health Expectancies ****************/ |
/*********** Health Expectancies ****************/ |
|
Line 5405 void concatwav(int wav[], int **dh, int
|
Line 5614 void concatwav(int wav[], int **dh, int
|
| xp[i] = x[i] + (i==theta ?delti[theta]:0); |
xp[i] = x[i] + (i==theta ?delti[theta]:0); |
| } |
} |
| |
|
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nresult); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 5437 void concatwav(int wav[], int **dh, int
|
Line 5646 void concatwav(int wav[], int **dh, int
|
| for(i=1; i<=npar; i++) /* Computes gradient x - delta */ |
for(i=1; i<=npar; i++) /* Computes gradient x - delta */ |
| xp[i] = x[i] - (i==theta ?delti[theta]:0); |
xp[i] = x[i] - (i==theta ?delti[theta]:0); |
| |
|
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nresult); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 5514 void concatwav(int wav[], int **dh, int
|
Line 5723 void concatwav(int wav[], int **dh, int
|
| /* end ppptj */ |
/* end ppptj */ |
| /* x centered again */ |
/* x centered again */ |
| |
|
| prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nresult); |
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nres); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 5609 void concatwav(int wav[], int **dh, int
|
Line 5818 void concatwav(int wav[], int **dh, int
|
| |
|
| pstamp(ficresvpl); |
pstamp(ficresvpl); |
| fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); |
fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); |
| fprintf(ficresvpl,"# Age"); |
fprintf(ficresvpl,"# Age "); |
| |
if(nresult >=1) |
| |
fprintf(ficresvpl," Result# "); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| fprintf(ficresvpl," %1d-%1d",i,i); |
fprintf(ficresvpl," %1d-%1d",i,i); |
| fprintf(ficresvpl,"\n"); |
fprintf(ficresvpl,"\n"); |
|
Line 5695 void concatwav(int wav[], int **dh, int
|
Line 5906 void concatwav(int wav[], int **dh, int
|
| varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ |
varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ |
| |
|
| fprintf(ficresvpl,"%.0f ",age ); |
fprintf(ficresvpl,"%.0f ",age ); |
| |
if(nresult >=1) |
| |
fprintf(ficresvpl,"%d ",nres ); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); |
fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); |
| fprintf(ficresvpl,"\n"); |
fprintf(ficresvpl,"\n"); |
|
Line 6099 void printinghtml(char fileresu[], char
|
Line 6312 void printinghtml(char fileresu[], char
|
| jj1=0; |
jj1=0; |
| |
|
| for(nres=1; nres <= nresult; nres++) /* For each resultline */ |
for(nres=1; nres <= nresult; nres++) /* For each resultline */ |
| for(k1=1; k1<=m;k1++){ |
for(k1=1; k1<=m;k1++){ /* For each combination of covariate */ |
| if(TKresult[nres]!= k1) |
if(TKresult[nres]!= k1) |
| continue; |
continue; |
| |
|
|
Line 6127 void printinghtml(char fileresu[], char
|
Line 6340 void printinghtml(char fileresu[], char
|
| } |
} |
| } |
} |
| /* aij, bij */ |
/* aij, bij */ |
| fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \ |
fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \ |
| <img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
<img src=\"%s_%d-1-%d.svg\">",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); |
| /* Pij */ |
/* Pij */ |
| fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \ |
fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2-%d.svg\">%s_%d-2-%d.svg</a><br> \ |
| <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
<img src=\"%s_%d-2-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); |
| /* Quasi-incidences */ |
/* Quasi-incidences */ |
| fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ |
fprintf(fichtm,"<br>\n- I<sub>ij</sub> 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, \ |
before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \ |
| incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \ |
incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \ |
| divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \ |
divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3-%d.svg\">%s_%d-3-%d.svg</a><br> \ |
| <img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); |
| /* Survival functions (period) in state j */ |
/* Survival functions (period) in state j */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1); |
<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); |
| } |
} |
| /* State specific survival functions (period) */ |
/* State specific survival functions (period) */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\ |
fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\ |
| Or probability to survive in various states (1 to %d) being in state %d at different ages. \ |
Or probability to survive in various states (1 to %d) being in state %d at different ages. \ |
| <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1); |
<a href=\"%s_%d-%d-%d.svg\">%s_%d%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); |
| } |
} |
| /* Period (stable) prevalence in each health state */ |
/* Period (stable) prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); |
<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); |
| } |
} |
| if(backcast==1){ |
if(backcast==1){ |
| /* Period (stable) back prevalence in each health state */ |
/* Period (stable) back prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1); |
<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); |
| } |
} |
| } |
} |
| if(prevfcast==1){ |
if(prevfcast==1){ |
| /* Projection of prevalence up to period (stable) prevalence in each health state */ |
/* Projection of prevalence up to period (stable) prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| <img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1); |
<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); |
| } |
} |
| } |
} |
| |
|
| for(cpt=1; cpt<=nlstate;cpt++) { |
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 in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \ |
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a> <br> \ |
| <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); |
<img src=\"%s_%d-%d-%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres); |
| } |
} |
| /* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
| }/* End k1 */ |
}/* End k1 */ |
|
Line 6230 See page 'Matrix of variance-covariance
|
Line 6443 See page 'Matrix of variance-covariance
|
| |
|
| jj1=0; |
jj1=0; |
| |
|
| for(nres=1; nres <= nresult; nres++) /* For each resultline */ |
for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
| for(k1=1; k1<=m;k1++){ |
for(k1=1; k1<=m;k1++){ |
| if(TKresult[nres]!= k1) |
if(TKresult[nres]!= k1) |
| continue; |
continue; |
|
Line 6254 See page 'Matrix of variance-covariance
|
Line 6467 See page 'Matrix of variance-covariance
|
| } |
} |
| for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \ |
fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \ |
| prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d.svg\"> %s_%d-%d.svg</a>\n <br>\ |
prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>\n <br>\ |
| <img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); |
<img src=\"%s_%d-%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
| } |
} |
| fprintf(fichtm,"\n<br>- Total life expectancy by age and \ |
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) \ |
health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ |
| true period expectancies (those weighted with period prevalences are also\ |
true period expectancies (those weighted with period prevalences are also\ |
| drawn in addition to the population based expectancies computed using\ |
drawn in addition to the population based expectancies computed using\ |
| observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\ |
observed and cahotic prevalences: <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>\n<br>\ |
| <img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); |
<img src=\"%s_%d-%d.svg\">",subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres); |
| /* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
| }/* End k1 */ |
}/* End k1 */ |
| |
}/* End nres */ |
| fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
| fflush(fichtm); |
fflush(fichtm); |
| } |
} |
|
Line 6330 void printinggnuplot(char fileresu[], ch
|
Line 6544 void printinggnuplot(char fileresu[], ch
|
| if(TKresult[nres]!= k1) |
if(TKresult[nres]!= k1) |
| continue; |
continue; |
| /* We are interested in selected combination by the resultline */ |
/* We are interested in selected combination by the resultline */ |
| printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); |
/* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */ |
| fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); |
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ |
|
Line 6339 void printinggnuplot(char fileresu[], ch
|
Line 6553 void printinggnuplot(char fileresu[], ch
|
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */ |
vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */ |
| /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */ |
/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */ |
| printf(" V%d=%d ",Tvaraff[k],vlv); |
/* printf(" V%d=%d ",Tvaraff[k],vlv); */ |
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ |
for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ |
| printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
/* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */ |
| fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
| } |
} |
| printf("\n#\n"); |
/* printf("\n#\n"); */ |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
| fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); |
fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); |
| fprintf(ficgp,"set xlabel \"Age\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); |
| set ylabel \"Probability\" \n \ |
|
| set ter svg size 640, 480\n \ |
|
| plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
|
| |
|
| for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); |
| for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); |
| for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
|
Line 6377 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 6588 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"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(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1)); |
| if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
| /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
| fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */ |
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */ |
| if(cptcoveff ==0){ |
if(cptcoveff ==0){ |
| fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ", 2+(cpt-1), cpt ); |
fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line lt 3", 2+(cpt-1), cpt ); |
| }else{ |
}else{ |
| kl=0; |
kl=0; |
| for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ |
for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ |
|
Line 6394 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 6605 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| /*6+1+(i-1)+(nlstate+1)*nlstate; 6+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*/ |
/* '' 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){ |
if(k==cptcoveff){ |
| 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], \ |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' w l lt 3",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ |
| 4+(cpt-1), cpt ); /* 4 or 6 ?*/ |
2+cptcoveff*2+(cpt-1), cpt ); /* 4 or 6 ?*/ |
| }else{ |
}else{ |
| fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]); |
fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]); |
| kl++; |
kl++; |
|
Line 6434 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 6645 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres); |
| for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| if(vpopbased==0) |
if(vpopbased==0) |
| fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); |
fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); |
|
Line 6464 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 6675 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); |
else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); |
| } /* state */ |
} /* state */ |
| } /* vpopbased */ |
} /* vpopbased */ |
| fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ |
fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */ |
| } /* end nres */ |
} /* end nres */ |
| } /* k1 end 2 eme*/ |
} /* k1 end 2 eme*/ |
| |
|
|
Line 6472 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 6683 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| /*3eme*/ |
/*3eme*/ |
| for (k1=1; k1<= m ; k1 ++){ |
for (k1=1; k1<= m ; k1 ++){ |
| for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
| if(TKresult[nres]!= k) |
if(TKresult[nres]!= k1) |
| continue; |
continue; |
| |
|
| for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
|
Line 6496 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 6707 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| |
|
| /* k=2+nlstate*(2*cpt-2); */ |
/* k=2+nlstate*(2*cpt-2); */ |
| k=2+(nlstate+1)*(cpt-1); |
k=2+(nlstate+1)*(cpt-1); |
| fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres); |
| fprintf(ficgp,"set ter svg size 640, 480\n\ |
fprintf(ficgp,"set ter svg size 640, 480\n\ |
| plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt); |
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"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); |
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); |
|
Line 6542 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 6753 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
| set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; |
k=3; |
|
Line 6588 set ter svg size 640, 480\nunset log y\n
|
Line 6799 set ter svg size 640, 480\nunset log y\n
|
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
| set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; |
k=3; |
|
Line 6643 set ter svg size 640, 480\nunset log y\n
|
Line 6854 set ter svg size 640, 480\nunset log y\n
|
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1,nres); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
| set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; /* Offset */ |
k=3; /* Offset */ |
|
Line 6689 set ter svg size 640, 480\nunset log y\n
|
Line 6900 set ter svg size 640, 480\nunset log y\n
|
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
| set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; /* Offset */ |
k=3; /* Offset */ |
|
Line 6741 set ter svg size 640, 480\nunset log y\n
|
Line 6952 set ter svg size 640, 480\nunset log y\n
|
| } |
} |
| |
|
| fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); |
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,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
| set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
| for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ |
for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ |
|
Line 6857 set ter svg size 640, 480\nunset log y\n
|
Line 7068 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng,nres); |
| fprintf(ficgp,"\nset ter svg size 640, 480 "); |
fprintf(ficgp,"\nset ter svg size 640, 480 "); |
| if (ng==1){ |
if (ng==1){ |
| fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ |
fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ |
|
Line 7832 int readdata(char datafile[], int firsto
|
Line 8043 int readdata(char datafile[], int firsto
|
| /*-------- data file ----------*/ |
/*-------- data file ----------*/ |
| FILE *fic; |
FILE *fic; |
| char dummy[]=" "; |
char dummy[]=" "; |
| int i=0, j=0, n=0, iv=0; |
int i=0, j=0, n=0, iv=0, v; |
| int lstra; |
int lstra; |
| int linei, month, year,iout; |
int linei, month, year,iout; |
| char line[MAXLINE], linetmp[MAXLINE]; |
char line[MAXLINE], linetmp[MAXLINE]; |
| char stra[MAXLINE], strb[MAXLINE]; |
char stra[MAXLINE], strb[MAXLINE]; |
| char *stratrunc; |
char *stratrunc; |
| |
|
| |
DummyV=ivector(1,NCOVMAX); /* 1 to 3 */ |
| |
FixedV=ivector(1,NCOVMAX); /* 1 to 3 */ |
| |
|
| |
for(v=1; v <=ncovcol;v++){ |
| |
DummyV[v]=0; |
| |
FixedV[v]=0; |
| |
} |
| |
for(v=ncovcol+1; v <=ncovcol+nqv;v++){ |
| |
DummyV[v]=1; |
| |
FixedV[v]=0; |
| |
} |
| |
for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){ |
| |
DummyV[v]=0; |
| |
FixedV[v]=1; |
| |
} |
| |
for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){ |
| |
DummyV[v]=1; |
| |
FixedV[v]=1; |
| |
} |
| |
for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){ |
| |
printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); |
| |
fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); |
| |
} |
| |
|
| if((fic=fopen(datafile,"r"))==NULL) { |
if((fic=fopen(datafile,"r"))==NULL) { |
| printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout); |
printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout); |
|
Line 8460 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 8693 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \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\ |
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); |
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;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} |
| for(v=1; v <=ncovcol;v++){ |
|
| DummyV[v]=0; |
|
| FixedV[v]=0; |
|
| } |
|
| for(v=ncovcol+1; v <=ncovcol+nqv;v++){ |
|
| DummyV[v]=1; |
|
| FixedV[v]=0; |
|
| } |
|
| for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){ |
|
| DummyV[v]=0; |
|
| FixedV[v]=1; |
|
| } |
|
| for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){ |
|
| DummyV[v]=1; |
|
| FixedV[v]=1; |
|
| } |
|
| for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){ |
|
| printf("Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); |
|
| fprintf(ficlog,"Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]); |
|
| } |
|
| for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ |
for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ |
| if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ |
if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ |
| Fixed[k]= 0; |
Fixed[k]= 0; |
|
Line 8505 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 8718 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| TvarFind[ncovf]=k; |
TvarFind[ncovf]=k; |
| TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
| TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
| }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */ |
}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */ |
| Fixed[k]= 0; |
Fixed[k]= 0; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| nqfveff++; |
nqfveff++; |
|
Line 8519 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 8732 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| TvarFind[ncovf]=k; |
TvarFind[ncovf]=k; |
| TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ |
TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ |
| TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ |
TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ |
| }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying variables */ |
}else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 0; |
Dummy[k]= 0; |
| ntveff++; /* Only simple time varying dummy variable */ |
ntveff++; /* Only simple time varying dummy variable */ |
|
Line 8530 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 8743 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| TvarsDind[nsd]=k; |
TvarsDind[nsd]=k; |
| ncovv++; /* Only simple time varying variables */ |
ncovv++; /* Only simple time varying variables */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; /* TvarVind[2]=2 TvarVind[3]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ |
| TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ |
TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ |
| TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ |
TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */ |
| printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv); |
printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv); |
|
Line 8546 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 8759 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| TvarsQ[nsq]=Tvar[k]; |
TvarsQ[nsq]=Tvar[k]; |
| TvarsQind[nsq]=k; |
TvarsQind[nsq]=k; |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ |
| TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ |
TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ |
| TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ |
TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ |
| TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */ |
TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */ |
|
Line 8558 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 8771 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| TvarA[ncova]=Tvar[k]; |
TvarA[ncova]=Tvar[k]; |
| TvarAind[ncova]=k; |
TvarAind[ncova]=k; |
| if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */ |
if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */ |
| Fixed[k]= 2; |
Fixed[k]= 2; |
| Dummy[k]= 2; |
Dummy[k]= 2; |
| modell[k].maintype= ATYPE; |
modell[k].maintype= ATYPE; |
| modell[k].subtype= APFD; |
modell[k].subtype= APFD; |
| /* ncoveff++; */ |
/* ncoveff++; */ |
| }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/ |
}else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/ |
| Fixed[k]= 2; |
Fixed[k]= 2; |
| Dummy[k]= 3; |
Dummy[k]= 3; |
| modell[k].maintype= ATYPE; |
modell[k].maintype= ATYPE; |
| modell[k].subtype= APFQ; /* Product age * fixed quantitative */ |
modell[k].subtype= APFQ; /* Product age * fixed quantitative */ |
| /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */ |
/* nqfveff++; /\* Only simple fixed quantitative variable *\/ */ |
| }else if( Tvar[k] <=ncovcol+nqv+ntv ){ |
}else if( Tvar[k] <=ncovcol+nqv+ntv ){ |
| Fixed[k]= 3; |
Fixed[k]= 3; |
| Dummy[k]= 2; |
Dummy[k]= 2; |
| modell[k].maintype= ATYPE; |
modell[k].maintype= ATYPE; |
| modell[k].subtype= APVD; /* Product age * varying dummy */ |
modell[k].subtype= APVD; /* Product age * varying dummy */ |
| /* ntveff++; /\* Only simple time varying dummy variable *\/ */ |
/* ntveff++; /\* Only simple time varying dummy variable *\/ */ |
| }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){ |
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){ |
| Fixed[k]= 3; |
Fixed[k]= 3; |
| Dummy[k]= 3; |
Dummy[k]= 3; |
| modell[k].maintype= ATYPE; |
modell[k].maintype= ATYPE; |
| modell[k].subtype= APVQ; /* Product age * varying quantitative */ |
modell[k].subtype= APVQ; /* Product age * varying quantitative */ |
| /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */ |
/* nqtveff++;/\* Only simple time varying quantitative variable *\/ */ |
| } |
} |
| }else if (Typevar[k] == 2) { /* product without age */ |
}else if (Typevar[k] == 2) { /* product without age */ |
| k1=Tposprod[k]; |
k1=Tposprod[k]; |
| if(Tvard[k1][1] <=ncovcol){ |
if(Tvard[k1][1] <=ncovcol){ |
| if(Tvard[k1][2] <=ncovcol){ |
if(Tvard[k1][2] <=ncovcol){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 0; |
Dummy[k]= 0; |
| modell[k].maintype= FTYPE; |
modell[k].maintype= FTYPE; |
| modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */ |
modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */ |
| ncovf++; /* Fixed variables without age */ |
ncovf++; /* Fixed variables without age */ |
| TvarF[ncovf]=Tvar[k]; |
TvarF[ncovf]=Tvar[k]; |
| TvarFind[ncovf]=k; |
TvarFind[ncovf]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv){ |
| Fixed[k]= 0; /* or 2 ?*/ |
Fixed[k]= 0; /* or 2 ?*/ |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= FTYPE; |
modell[k].maintype= FTYPE; |
| modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */ |
modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */ |
| ncovf++; /* Varying variables without age */ |
ncovf++; /* Varying variables without age */ |
| TvarF[ncovf]=Tvar[k]; |
TvarF[ncovf]=Tvar[k]; |
| TvarFind[ncovf]=k; |
TvarFind[ncovf]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 0; |
Dummy[k]= 0; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */ |
modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */ |
modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| } |
} |
| }else if(Tvard[k1][1] <=ncovcol+nqv){ |
}else if(Tvard[k1][1] <=ncovcol+nqv){ |
| if(Tvard[k1][2] <=ncovcol){ |
if(Tvard[k1][2] <=ncovcol){ |
| Fixed[k]= 0; /* or 2 ?*/ |
Fixed[k]= 0; /* or 2 ?*/ |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= FTYPE; |
modell[k].maintype= FTYPE; |
| modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */ |
modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */ |
| ncovf++; /* Fixed variables without age */ |
ncovf++; /* Fixed variables without age */ |
| TvarF[ncovf]=Tvar[k]; |
TvarF[ncovf]=Tvar[k]; |
| TvarFind[ncovf]=k; |
TvarFind[ncovf]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */ |
modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */ |
modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| } |
} |
| }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ |
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ |
| if(Tvard[k1][2] <=ncovcol){ |
if(Tvard[k1][2] <=ncovcol){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */ |
modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */ |
modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 0; |
Dummy[k]= 0; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */ |
modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */ |
modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| } |
} |
| }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ |
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ |
| if(Tvard[k1][2] <=ncovcol){ |
if(Tvard[k1][2] <=ncovcol){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */ |
modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */ |
modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */ |
modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
| Fixed[k]= 1; |
Fixed[k]= 1; |
| Dummy[k]= 1; |
Dummy[k]= 1; |
| modell[k].maintype= VTYPE; |
modell[k].maintype= VTYPE; |
| modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */ |
modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */ |
| ncovv++; /* Varying variables without age */ |
ncovv++; /* Varying variables without age */ |
| TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
| TvarVind[ncovv]=k; |
TvarVind[ncovv]=k; |
| } |
} |
| }else{ |
}else{ |
| printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); |
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]); |
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 */ |
} /*end k1*/ |
| }else{ |
}else{ |
| printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); |
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]); |
fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); |
|
Line 8760 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 8973 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| } |
} |
| |
|
| int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) |
int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) |
| { |
{/* Check ages at death */ |
| int i, m; |
int i, m; |
| int firstone=0; |
int firstone=0; |
| |
|
|
Line 9233 int back_prevalence_limit(double *p, dou
|
Line 9446 int back_prevalence_limit(double *p, dou
|
| if(mobilavproj > 0){ |
if(mobilavproj > 0){ |
| /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */ |
/* 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, 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, nres); |
| }else if (mobilavproj == 0){ |
}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); |
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); |
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); |
exit(1); |
| }else{ |
}else{ |
| /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, 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); |
bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k,nres); |
| } |
} |
| fprintf(ficresplb,"%.0f ",age ); |
fprintf(ficresplb,"%.0f ",age ); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
|
Line 9493 int main(int argc, char *argv[])
|
Line 9706 int main(int argc, char *argv[])
|
| double **prlim; |
double **prlim; |
| double **bprlim; |
double **bprlim; |
| double ***param; /* Matrix of parameters */ |
double ***param; /* Matrix of parameters */ |
| double *p; |
double ***paramstart; /* Matrix of starting parameter values */ |
| |
double *p, *pstart; /* p=param[1][1] pstart is for starting values guessed by freqsummary */ |
| double **matcov; /* Matrix of covariance */ |
double **matcov; /* Matrix of covariance */ |
| double **hess; /* Hessian matrix */ |
double **hess; /* Hessian matrix */ |
| double ***delti3; /* Scale */ |
double ***delti3; /* Scale */ |
|
Line 9804 int main(int argc, char *argv[])
|
Line 10018 int main(int argc, char *argv[])
|
| delti=delti3[1][1]; |
delti=delti3[1][1]; |
| /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ |
/*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ |
| if(mle==-1){ /* Print a wizard for help writing covariance matrix */ |
if(mle==-1){ /* Print a wizard for help writing covariance matrix */ |
| |
/* We could also provide initial parameters values giving by simple logistic regression |
| |
* only one way, that is without matrix product. We will have nlstate maximizations */ |
| |
/* for(i=1;i<nlstate;i++){ */ |
| |
/* /\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */ |
| |
/* mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */ |
| |
/* } */ |
| prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
| printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
| fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
|
Line 9833 int main(int argc, char *argv[])
|
Line 10053 int main(int argc, char *argv[])
|
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
| |
paramstart= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
| for(i=1; i <=nlstate; i++){ |
for(i=1; i <=nlstate; i++){ |
| j=0; |
j=0; |
| for(jj=1; jj <=nlstate+ndeath; jj++){ |
for(jj=1; jj <=nlstate+ndeath; jj++){ |
|
Line 9869 run imach with mle=-1 to get a correct t
|
Line 10090 run imach with mle=-1 to get a correct t
|
| } |
} |
| fflush(ficlog); |
fflush(ficlog); |
| |
|
| /* Reads scales values */ |
/* Reads parameters values */ |
| p=param[1][1]; |
p=param[1][1]; |
| |
pstart=paramstart[1][1]; |
| |
|
| /* Reads comments: lines beginning with '#' */ |
/* Reads comments: lines beginning with '#' */ |
| while((c=getc(ficpar))=='#' && c!= EOF){ |
while((c=getc(ficpar))=='#' && c!= EOF){ |
|
Line 10050 Please run with mle=-1 to get a correct
|
Line 10272 Please run with mle=-1 to get a correct
|
| Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */ |
Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */ |
| Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */ |
Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */ |
| Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */ |
Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */ |
| DummyV=ivector(1,NCOVMAX); /* 1 to 3 */ |
|
| FixedV=ivector(1,NCOVMAX); /* 1 to 3 */ |
|
| /* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs). |
/* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs). |
| For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4, |
For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4, |
| Tvar[4=age*V3] is 3 and 'age' is recorded in Tage. |
Tvar[4=age*V3] is 3 and 'age' is recorded in Tage. |
|
Line 10137 Please run with mle=-1 to get a correct
|
Line 10357 Please run with mle=-1 to get a correct
|
| */ |
*/ |
| |
|
| concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); |
concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); |
| /* */ |
/* Concatenates waves */ |
| |
|
| free_vector(moisdc,1,n); |
free_vector(moisdc,1,n); |
| free_vector(andc,1,n); |
free_vector(andc,1,n); |
|
Line 10318 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 10538 Title=%s <br>Datafile=%s Firstpass=%d La
|
| /* Calculates basic frequencies. Computes observed prevalence at single age |
/* Calculates basic frequencies. Computes observed prevalence at single age |
| and for any valid combination of covariates |
and for any valid combination of covariates |
| and prints on file fileres'p'. */ |
and prints on file fileres'p'. */ |
| freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ |
freqsummary(fileres, p, pstart, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ |
| firstpass, lastpass, stepm, weightopt, model); |
firstpass, lastpass, stepm, weightopt, model); |
| |
|
| fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
|
Line 10340 Interval (in months) between two waves:
|
Line 10560 Interval (in months) between two waves:
|
| /* For mortality only */ |
/* For mortality only */ |
| if (mle==-3){ |
if (mle==-3){ |
| ximort=matrix(1,NDIM,1,NDIM); |
ximort=matrix(1,NDIM,1,NDIM); |
| for(i=1;i<=NDIM;i++) |
for(i=1;i<=NDIM;i++) |
| for(j=1;j<=NDIM;j++) |
for(j=1;j<=NDIM;j++) |
| ximort[i][j]=0.; |
ximort[i][j]=0.; |
| /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ |
/* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ |
| cens=ivector(1,n); |
cens=ivector(1,n); |
| ageexmed=vector(1,n); |
ageexmed=vector(1,n); |
|
Line 10578 Please run with mle=-1 to get a correct
|
Line 10798 Please run with mle=-1 to get a correct
|
| printf("\n"); |
printf("\n"); |
| if(mle>=1){ /* Could be 1 or 2, Real Maximization */ |
if(mle>=1){ /* Could be 1 or 2, Real Maximization */ |
| /* mlikeli uses func not funcone */ |
/* mlikeli uses func not funcone */ |
| |
/* for(i=1;i<nlstate;i++){ */ |
| |
/* /\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */ |
| |
/* mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */ |
| |
/* } */ |
| mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); |
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); |
| } |
} |
| if(mle==0) {/* No optimization, will print the likelihoods for the datafile */ |
if(mle==0) {/* No optimization, will print the likelihoods for the datafile */ |
|
Line 10869 Please run with mle=-1 to get a correct
|
Line 11093 Please run with mle=-1 to get a correct
|
| }else |
}else |
| break; |
break; |
| } |
} |
| |
if (!feof(ficpar)) |
| while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){ |
while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){ |
| if (num_filled == 0) |
if (num_filled == 0){ |
| resultline[0]='\0'; |
resultline[0]='\0'; |
| else if (num_filled != 1){ |
break; |
| |
} else if (num_filled != 1){ |
| printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line); |
printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line); |
| } |
} |
| nresult++; /* Sum of resultlines */ |
nresult++; /* Sum of resultlines */ |
|
Line 10902 Please run with mle=-1 to get a correct
|
Line 11128 Please run with mle=-1 to get a correct
|
| break; |
break; |
| else{ /* Processess output results for this combination of covariate values */ |
else{ /* Processess output results for this combination of covariate values */ |
| } |
} |
| } |
} /* end while */ |
| |
|
| |
|
| |
|
|
Line 11120 Please run with mle=-1 to get a correct
|
Line 11346 Please run with mle=-1 to get a correct
|
| for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ |
for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ |
| if(TKresult[nres]!= k) |
if(TKresult[nres]!= k) |
| continue; |
continue; |
| printf("\n#****** Selected:"); |
printf("\n#****** Result for:"); |
| fprintf(ficrest,"\n#****** Selected:"); |
fprintf(ficrest,"\n#****** Result for:"); |
| fprintf(ficlog,"\n#****** Selected:"); |
fprintf(ficlog,"\n#****** Result for:"); |
| for(j=1;j<=cptcoveff;j++){ |
for(j=1;j<=cptcoveff;j++){ |
| printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,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(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src