|
version 1.203, 2015/09/30 17:45:14
|
version 1.215, 2015/12/16 08:52:24
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.215 2015/12/16 08:52:24 brouard |
| |
Summary: 0.98r4 working |
| |
|
| |
Revision 1.214 2015/12/16 06:57:54 brouard |
| |
Summary: temporary not working |
| |
|
| |
Revision 1.213 2015/12/11 18:22:17 brouard |
| |
Summary: 0.98r4 |
| |
|
| |
Revision 1.212 2015/11/21 12:47:24 brouard |
| |
Summary: minor typo |
| |
|
| |
Revision 1.211 2015/11/21 12:41:11 brouard |
| |
Summary: 0.98r3 with some graph of projected cross-sectional |
| |
|
| |
Author: Nicolas Brouard |
| |
|
| |
Revision 1.210 2015/11/18 17:41:20 brouard |
| |
Summary: Start working on projected prevalences |
| |
|
| |
Revision 1.209 2015/11/17 22:12:03 brouard |
| |
Summary: Adding ftolpl parameter |
| |
Author: N Brouard |
| |
|
| |
We had difficulties to get smoothed confidence intervals. It was due |
| |
to the period prevalence which wasn't computed accurately. The inner |
| |
parameter ftolpl is now an outer parameter of the .imach parameter |
| |
file after estepm. If ftolpl is small 1.e-4 and estepm too, |
| |
computation are long. |
| |
|
| |
Revision 1.208 2015/11/17 14:31:57 brouard |
| |
Summary: temporary |
| |
|
| |
Revision 1.207 2015/10/27 17:36:57 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.206 2015/10/24 07:14:11 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.205 2015/10/23 15:50:53 brouard |
| |
Summary: 0.98r3 some clarification for graphs on likelihood contributions |
| |
|
| |
Revision 1.204 2015/10/01 16:20:26 brouard |
| |
Summary: Some new graphs of contribution to likelihood |
| |
|
| Revision 1.203 2015/09/30 17:45:14 brouard |
Revision 1.203 2015/09/30 17:45:14 brouard |
| Summary: looking at better estimation of the hessian |
Summary: looking at better estimation of the hessian |
| |
|
|
Line 729 typedef struct {
|
Line 774 typedef struct {
|
| #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
| #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
| #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| |
/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/ |
| |
#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 |
| #define MAXN 20000 |
#define MAXN 20000 |
| #define YEARM 12. /**< Number of months per year */ |
#define YEARM 12. /**< Number of months per year */ |
| #define AGESUP 130 |
#define AGESUP 130 |
|
Line 749 typedef struct {
|
Line 796 typedef struct {
|
| /* $State$ */ |
/* $State$ */ |
| #include "version.h" |
#include "version.h" |
| char version[]=__IMACH_VERSION__; |
char version[]=__IMACH_VERSION__; |
| char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; |
char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; |
| char fullversion[]="$Revision$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
| char strstart[80]; |
char strstart[80]; |
| char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
|
Line 786 double **matprod2(); /* test */
|
Line 833 double **matprod2(); /* test */
|
| double **oldm, **newm, **savm; /* Working pointers to matrices */ |
double **oldm, **newm, **savm; /* Working pointers to matrices */ |
| double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ |
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ |
| /*FILE *fic ; */ /* Used in readdata only */ |
/*FILE *fic ; */ /* Used in readdata only */ |
| FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; |
FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; |
| FILE *ficlog, *ficrespow; |
FILE *ficlog, *ficrespow; |
| int globpr=0; /* Global variable for printing or not */ |
int globpr=0; /* Global variable for printing or not */ |
| double fretone; /* Only one call to likelihood */ |
double fretone; /* Only one call to likelihood */ |
|
Line 940 static int split( char *path, char *dirc
|
Line 987 static int split( char *path, char *dirc
|
| } |
} |
| /* got dirc from getcwd*/ |
/* got dirc from getcwd*/ |
| printf(" DIRC = %s \n",dirc); |
printf(" DIRC = %s \n",dirc); |
| } else { /* strip direcotry from path */ |
} else { /* strip directory from path */ |
| ss++; /* after this, the filename */ |
ss++; /* after this, the filename */ |
| l2 = strlen( ss ); /* length of filename */ |
l2 = strlen( ss ); /* length of filename */ |
| if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); |
if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); |
|
Line 1348 char *subdirf3(char fileres[], char *pre
|
Line 1395 char *subdirf3(char fileres[], char *pre
|
| strcat(tmpout,fileres); |
strcat(tmpout,fileres); |
| return tmpout; |
return tmpout; |
| } |
} |
| |
|
| |
/*************** function subdirfext ***********/ |
| |
char *subdirfext(char fileres[], char *preop, char *postop) |
| |
{ |
| |
|
| |
strcpy(tmpout,preop); |
| |
strcat(tmpout,fileres); |
| |
strcat(tmpout,postop); |
| |
return tmpout; |
| |
} |
| |
|
| |
/*************** function subdirfext3 ***********/ |
| |
char *subdirfext3(char fileres[], char *preop, char *postop) |
| |
{ |
| |
|
| |
/* Caution optionfilefiname is hidden */ |
| |
strcpy(tmpout,optionfilefiname); |
| |
strcat(tmpout,"/"); |
| |
strcat(tmpout,preop); |
| |
strcat(tmpout,fileres); |
| |
strcat(tmpout,postop); |
| |
return tmpout; |
| |
} |
| |
|
| char *asc_diff_time(long time_sec, char ascdiff[]) |
char *asc_diff_time(long time_sec, char ascdiff[]) |
| { |
{ |
| long sec_left, days, hours, minutes; |
long sec_left, days, hours, minutes; |
|
Line 1954 double **prevalim(double **prlim, int nl
|
Line 2024 double **prevalim(double **prlim, int nl
|
| { |
{ |
| /* Computes the prevalence limit in each live state at age x by left multiplying the unit |
/* Computes the prevalence limit in each live state at age x 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 */ |
| |
/* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ |
| |
/* Wx is row vector: population in state 1, population in state 2, population dead */ |
| |
/* or prevalence in state 1, prevalence in state 2, 0 */ |
| |
/* newm is the matrix after multiplications, its rows are identical at a factor */ |
| |
/* Initial matrix pimij */ |
| |
/* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ |
| |
/* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ |
| |
/* 0, 0 , 1} */ |
| |
/* |
| |
* and after some iteration: */ |
| |
/* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */ |
| |
/* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */ |
| |
/* 0, 0 , 1} */ |
| |
/* And prevalence by suppressing the deaths are close to identical rows in prlim: */ |
| |
/* {0.51571254859325999, 0.4842874514067399, */ |
| |
/* 0.51326036147820708, 0.48673963852179264} */ |
| |
/* If we start from prlim again, prlim tends to a constant matrix */ |
| |
|
| int i, ii,j,k; |
int i, ii,j,k; |
| double min, max, maxmin, maxmax,sumnew=0.; |
double *min, *max, *meandiff, maxmax,sumnew=0.; |
| /* double **matprod2(); */ /* test */ |
/* double **matprod2(); */ /* test */ |
| double **out, cov[NCOVMAX+1], **pmij(); |
double **out, cov[NCOVMAX+1], **pmij(); |
| double **newm; |
double **newm; |
| double agefin, delaymax=100 ; /* Max number of years to converge */ |
double agefin, delaymax=200. ; /* 100 Max number of years to converge */ |
| int ncvloop=0; |
int ncvloop=0; |
| |
|
| |
min=vector(1,nlstate); |
| |
max=vector(1,nlstate); |
| |
meandiff=vector(1,nlstate); |
| |
|
| 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); |
|
Line 2000 double **prevalim(double **prlim, int nl
|
Line 2091 double **prevalim(double **prlim, int nl
|
| |
|
| savm=oldm; |
savm=oldm; |
| oldm=newm; |
oldm=newm; |
| maxmax=0.; |
|
| for(j=1;j<=nlstate;j++){ |
for(j=1; j<=nlstate; j++){ |
| min=1.; |
max[j]=0.; |
| max=0.; |
min[j]=1.; |
| for(i=1; i<=nlstate; i++) { |
} |
| sumnew=0; |
for(i=1;i<=nlstate;i++){ |
| for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; |
sumnew=0; |
| |
for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; |
| |
for(j=1; j<=nlstate; j++){ |
| prlim[i][j]= newm[i][j]/(1-sumnew); |
prlim[i][j]= newm[i][j]/(1-sumnew); |
| max=FMAX(max,prlim[i][j]); |
max[j]=FMAX(max[j],prlim[i][j]); |
| min=FMIN(min,prlim[i][j]); |
min[j]=FMIN(min[j],prlim[i][j]); |
| /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */ |
|
| } |
} |
| maxmin=(max-min)/(max+min)*2; |
} |
| maxmax=FMAX(maxmax,maxmin); |
|
| |
maxmax=0.; |
| |
for(j=1; j<=nlstate; j++){ |
| |
meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */ |
| |
maxmax=FMAX(maxmax,meandiff[j]); |
| |
/* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */ |
| } /* j loop */ |
} /* j loop */ |
| *ncvyear= (int)age- (int)agefin; |
*ncvyear= (int)age- (int)agefin; |
| /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
/* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
| if(maxmax < ftolpl){ |
if(maxmax < ftolpl){ |
| /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
/* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
| |
free_vector(min,1,nlstate); |
| |
free_vector(max,1,nlstate); |
| |
free_vector(meandiff,1,nlstate); |
| return prlim; |
return prlim; |
| } |
} |
| } /* age loop */ |
} /* age loop */ |
| printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g. \n\ |
/* After some age loop it doesn't converge */ |
| Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ |
| /* 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); */ |
Earliest 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); */ |
| |
free_vector(min,1,nlstate); |
| |
free_vector(max,1,nlstate); |
| |
free_vector(meandiff,1,nlstate); |
| |
|
| return prlim; /* should not reach here */ |
return prlim; /* should not reach here */ |
| } |
} |
| |
|
|
Line 2150 double ***hpxij(double ***po, int nhstep
|
Line 2255 double ***hpxij(double ***po, int nhstep
|
| double **out, cov[NCOVMAX+1]; |
double **out, cov[NCOVMAX+1]; |
| double **newm; |
double **newm; |
| double agexact; |
double agexact; |
| |
double agebegin, ageend; |
| |
|
| /* Hstepm could be zero and should return the unit matrix */ |
/* Hstepm could be zero and should return the unit matrix */ |
| for (i=1;i<=nlstate+ndeath;i++) |
for (i=1;i<=nlstate+ndeath;i++) |
|
Line 2163 double ***hpxij(double ***po, int nhstep
|
Line 2269 double ***hpxij(double ***po, int nhstep
|
| newm=savm; |
newm=savm; |
| /* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
| cov[1]=1.; |
cov[1]=1.; |
| agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; |
agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ |
| cov[2]=agexact; |
cov[2]=agexact; |
| if(nagesqr==1) |
if(nagesqr==1) |
| cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
|
Line 2534 double funcone( double *x)
|
Line 2640 double funcone( double *x)
|
| int s1, s2; |
int s1, s2; |
| double bbh, survp; |
double bbh, survp; |
| double agexact; |
double agexact; |
| |
double agebegin, ageend; |
| /*extern weight */ |
/*extern weight */ |
| /* We are differentiating ll according to initial status */ |
/* We are differentiating ll according to initial status */ |
| /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
|
Line 2552 double funcone( double *x)
|
Line 2659 double funcone( double *x)
|
| 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); |
| } |
} |
| for(d=0; d<dh[mi][i]; d++){ |
|
| |
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 */ |
| |
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] |
| |
and mw[mi+1][i]. dh depends on stepm.*/ |
| newm=savm; |
newm=savm; |
| agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
| cov[2]=agexact; |
cov[2]=agexact; |
|
Line 2600 double funcone( double *x)
|
Line 2712 double funcone( double *x)
|
| 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 %6d %2d %2d %2d %2d %3d %11.6f %8.4f\ |
fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ |
| %11.6f %11.6f %11.6f ", \ |
%11.6f %11.6f %11.6f ", \ |
| num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], |
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; |
|
Line 2640 void likelione(FILE *ficres,double p[],
|
Line 2752 void likelione(FILE *ficres,double p[],
|
| printf("Problem with resultfile: %s\n", fileresilk); |
printf("Problem with resultfile: %s\n", fileresilk); |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); |
| } |
} |
| fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); |
fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); |
| fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight 2wlli out sav "); |
fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav "); |
| /* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ |
/* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ |
| for(k=1; k<=nlstate; k++) |
for(k=1; k<=nlstate; k++) |
| fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); |
fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); |
|
Line 2651 void likelione(FILE *ficres,double p[],
|
Line 2763 void likelione(FILE *ficres,double p[],
|
| *fretone=(*funcone)(p); |
*fretone=(*funcone)(p); |
| if(*globpri !=0){ |
if(*globpri !=0){ |
| fclose(ficresilk); |
fclose(ficresilk); |
| fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with initial parameters and mle >= 1. You should at least run with mle >= 1 and starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk)); |
if (mle ==0) |
| fprintf(fichtm,"<br>- The first 3 individuals are drawn with lines. The function drawn is -2Log(L) in log scale: <a href=\"%s.png\">%s.png</a><br> \ |
fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with initial parameters and mle = %d.",mle); |
| <img src=\"%s.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); |
else if(mle >=1) |
| |
fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); |
| |
fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk)); |
| |
|
| |
|
| |
for (k=1; k<= nlstate ; k++) { |
| |
fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ |
| |
<img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); |
| |
} |
| |
fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \ |
| |
<img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); |
| |
fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \ |
| |
<img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); |
| fflush(fichtm); |
fflush(fichtm); |
| } |
} |
| return; |
return; |
| } |
} |
| |
|
|
Line 2919 double hessij( double x[], double **hess
|
Line 3043 double hessij( double x[], double **hess
|
| double p2[MAXPARM+1]; |
double p2[MAXPARM+1]; |
| int k, kmax=1; |
int k, kmax=1; |
| double v1, v2, cv12, lc1, lc2; |
double v1, v2, cv12, lc1, lc2; |
| |
|
| |
int firstime=0; |
| |
|
| fx=func(x); |
fx=func(x); |
| for (k=1; k<=kmax; k=k+10) { |
for (k=1; k<=kmax; k=k+10) { |
|
Line 2940 double hessij( double x[], double **hess
|
Line 3066 double hessij( double x[], double **hess
|
| k4=func(p2)-fx; |
k4=func(p2)-fx; |
| res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */ |
res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */ |
| if(k1*k2*k3*k4 <0.){ |
if(k1*k2*k3*k4 <0.){ |
| |
firstime=1; |
| kmax=kmax+10; |
kmax=kmax+10; |
| if(kmax >=10){ |
} |
| |
if(kmax >=10 || firstime ==1){ |
| printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); |
printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); |
| fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); |
fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; 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); |
| } |
|
| } |
} |
| #ifdef DEBUGHESSIJ |
#ifdef DEBUGHESSIJ |
| v1=hess[thetai][thetai]; |
v1=hess[thetai][thetai]; |
|
Line 3100 void pstamp(FILE *fichier)
|
Line 3227 void pstamp(FILE *fichier)
|
| } |
} |
| |
|
| /************ Frequencies ********************/ |
/************ Frequencies ********************/ |
| void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) |
void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ |
| |
int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],\ |
| |
int firstpass, int lastpass, int stepm, int weightopt, char model[]) |
| { /* Some frequencies */ |
{ /* Some frequencies */ |
| |
|
| int i, m, jk, j1, bool, z1,j; |
int i, m, jk, j1, bool, z1,j; |
| |
int mi; /* Effective wave */ |
| int first; |
int first; |
| double ***freq; /* Frequencies */ |
double ***freq; /* Frequencies */ |
| double *pp, **prop; |
double *pp, **prop; |
| double pos,posprop, k2, dateintsum=0,k2cpt=0; |
double pos,posprop, k2, dateintsum=0,k2cpt=0; |
| char fileresp[FILENAMELENGTH]; |
char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; |
| |
double agebegin, ageend; |
| |
|
| pp=vector(1,nlstate); |
pp=vector(1,nlstate); |
| prop=matrix(1,nlstate,iagemin,iagemax+3); |
prop=matrix(1,nlstate,iagemin,iagemax+3); |
| strcpy(fileresp,"P_"); |
strcpy(fileresp,"P_"); |
| strcat(fileresp,fileresu); |
strcat(fileresp,fileresu); |
| |
/*strcat(fileresphtm,fileresu);*/ |
| if((ficresp=fopen(fileresp,"w"))==NULL) { |
if((ficresp=fopen(fileresp,"w"))==NULL) { |
| printf("Problem with prevalence resultfile: %s\n", fileresp); |
printf("Problem with prevalence resultfile: %s\n", fileresp); |
| 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")); |
| |
if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { |
| |
printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
| |
fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
| |
fflush(ficlog); |
| |
exit(70); |
| |
} |
| |
else{ |
| |
fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
| |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
| |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
| |
fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
| |
} |
| |
fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm); |
| |
|
| |
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); |
| |
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { |
| |
printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
| |
fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
| |
fflush(ficlog); |
| |
exit(70); |
| |
} |
| |
else{ |
| |
fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
| |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
| |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
| |
fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
| |
} |
| |
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr); |
| |
|
| freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); |
freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); |
| j1=0; |
j1=0; |
| |
|
|
Line 3127 void freqsummary(char fileres[], int ia
|
Line 3290 void freqsummary(char fileres[], int ia
|
| |
|
| first=1; |
first=1; |
| |
|
| /* for(k1=1; k1<=j ; k1++){ */ /* Loop on covariates */ |
for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */ |
| /* for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */ |
|
| /* j1++; */ |
|
| for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ |
|
| /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
| scanf("%d", i);*/ |
scanf("%d", i);*/ |
| for (i=-5; i<=nlstate+ndeath; i++) |
for (i=-5; i<=nlstate+ndeath; i++) |
|
Line 3144 void freqsummary(char fileres[], int ia
|
Line 3304 void freqsummary(char fileres[], int ia
|
| |
|
| dateintsum=0; |
dateintsum=0; |
| k2cpt=0; |
k2cpt=0; |
| for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { /* For each individual i */ |
| bool=1; |
bool=1; |
| if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
| for (z1=1; z1<=cptcoveff; z1++) |
for (z1=1; z1<=cptcoveff; z1++) |
| if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
| /* Tests if the value of each of the covariates of i is equal to filter j1 */ |
/* Tests if the value of each of the covariates of i is equal to filter j1 */ |
|
Line 3156 void freqsummary(char fileres[], int ia
|
Line 3316 void freqsummary(char fileres[], int ia
|
| j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
| /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
| } |
} |
| } |
} /* cptcovn > 0 */ |
| |
|
| if (bool==1){ |
if (bool==1){ |
| for(m=firstpass; m<=lastpass; m++){ |
/* for(m=firstpass; m<=lastpass; m++){ */ |
| k2=anint[m][i]+(mint[m][i]/12.); |
for(mi=1; mi<wav[i];mi++){ |
| /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ |
m=mw[mi][i]; |
| if(agev[m][i]==0) agev[m][i]=iagemax+1; |
/* dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective (mi) waves m=mw[mi][i] |
| if(agev[m][i]==1) agev[m][i]=iagemax+2; |
and mw[mi+1][i]. dh depends on stepm. */ |
| if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; |
agebegin=agev[m][i]; /* Age at beginning of wave before transition*/ |
| |
ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /* Age at end of wave and transition */ |
| |
if(m >=firstpass && m <=lastpass){ |
| |
k2=anint[m][i]+(mint[m][i]/12.); |
| |
/*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ |
| |
if(agev[m][i]==0) agev[m][i]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ |
| |
if(agev[m][i]==1) agev[m][i]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ |
| |
if (s[m][i]>0 && s[m][i]<=nlstate) /* If status at wave m is known and a live state */ |
| |
prop[s[m][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */ |
| if (m<lastpass) { |
if (m<lastpass) { |
| freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; |
/* if(s[m][i]==4 && s[m+1][i]==4) */ |
| freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; |
/* printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i]); */ |
| } |
if(s[m][i]==-1) |
| |
printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i],agebegin, ageend, (int)((agebegin+ageend)/2.)); |
| if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3))) { |
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */ |
| dateintsum=dateintsum+k2; |
/* freq[s[m][i]][s[m+1][i]][(int)((agebegin+ageend)/2.)] += weight[i]; */ |
| k2cpt++; |
freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */ |
| } |
} |
| /*}*/ |
} |
| } |
if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) { |
| } |
dateintsum=dateintsum+k2; |
| } /* end i */ |
k2cpt++; |
| |
/* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ |
| |
} |
| |
/*}*/ |
| |
} /* end m */ |
| |
} /* end bool */ |
| |
} /* end i = 1 to imx */ |
| |
|
| /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
| pstamp(ficresp); |
pstamp(ficresp); |
| if (cptcovn>0) { |
if (cptcovn>0) { |
| fprintf(ficresp, "\n#********** Variable "); |
fprintf(ficresp, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
| fprintf(ficresp, "**********\n#"); |
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); |
| |
for (z1=1; z1<=cptcoveff; z1++){ |
| |
fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
} |
| |
fprintf(ficresp, "**********\n#"); |
| |
fprintf(ficresphtm, "**********</h3>\n"); |
| |
fprintf(ficresphtmfr, "**********</h3>\n"); |
| fprintf(ficlog, "\n#********** Variable "); |
fprintf(ficlog, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| fprintf(ficlog, "**********\n#"); |
fprintf(ficlog, "**********\n"); |
| } |
} |
| for(i=1; i<=nlstate;i++) |
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(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(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(i=iagemin; i <= iagemax+3; i++){ |
for(i=iagemin; i <= iagemax+3; i++){ |
| if(i==iagemax+3){ |
fprintf(ficresphtm,"<tr>"); |
| |
if(i==iagemax+1){ |
| |
fprintf(ficlog,"1"); |
| |
fprintf(ficresphtmfr,"<tr><th>0</th> "); |
| |
}else if(i==iagemax+2){ |
| |
fprintf(ficlog,"0"); |
| |
fprintf(ficresphtmfr,"<tr><th>Unknown</th> "); |
| |
}else if(i==iagemax+3){ |
| fprintf(ficlog,"Total"); |
fprintf(ficlog,"Total"); |
| |
fprintf(ficresphtmfr,"<tr><th>Total</th> "); |
| }else{ |
}else{ |
| if(first==1){ |
if(first==1){ |
| first=0; |
first=0; |
| printf("See log file for details...\n"); |
printf("See log file for details...\n"); |
| } |
} |
| |
fprintf(ficresphtmfr,"<tr><th>%d</th> ",i); |
| fprintf(ficlog,"Age %d", i); |
fprintf(ficlog,"Age %d", i); |
| } |
} |
| for(jk=1; jk <=nlstate ; jk++){ |
for(jk=1; jk <=nlstate ; jk++){ |
|
Line 3243 void freqsummary(char fileres[], int ia
|
Line 3450 void freqsummary(char fileres[], int ia
|
| if( i <= iagemax){ |
if( i <= iagemax){ |
| if(pos>=1.e-5){ |
if(pos>=1.e-5){ |
| fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); |
fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); |
| |
fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",i,prop[jk][i]/posprop, prop[jk][i],posprop); |
| /*probs[i][jk][j1]= pp[jk]/pos;*/ |
/*probs[i][jk][j1]= pp[jk]/pos;*/ |
| /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ |
/*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ |
| } |
} |
| else |
else{ |
| fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); |
fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); |
| |
fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",i, prop[jk][i],posprop); |
| |
} |
| } |
} |
| } |
} |
| |
|
| 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][i] !=0 ) { |
if(freq[jk][m][i] !=0 ) { /* minimizing output */ |
| if(first==1) |
if(first==1){ |
| printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); |
printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); |
| |
} |
| fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]); |
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]); |
| } |
} |
| if(i <= iagemax) |
if(jk!=0 && m!=0) |
| |
fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][i]); |
| |
} |
| |
} |
| |
fprintf(ficresphtmfr,"</tr>\n "); |
| |
if(i <= iagemax){ |
| fprintf(ficresp,"\n"); |
fprintf(ficresp,"\n"); |
| |
fprintf(ficresphtm,"</tr>\n"); |
| |
} |
| if(first==1) |
if(first==1) |
| printf("Others in log...\n"); |
printf("Others in log...\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| } |
} /* end loop i */ |
| |
fprintf(ficresphtm,"</table>\n"); |
| |
fprintf(ficresphtmfr,"</table>\n"); |
| /*}*/ |
/*}*/ |
| } |
} /* end j1 */ |
| dateintmean=dateintsum/k2cpt; |
dateintmean=dateintsum/k2cpt; |
| |
|
| fclose(ficresp); |
fclose(ficresp); |
| |
fclose(ficresphtm); |
| |
fclose(ficresphtmfr); |
| free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); |
free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); |
| free_vector(pp,1,nlstate); |
free_vector(pp,1,nlstate); |
| free_matrix(prop,1,nlstate,iagemin, iagemax+3); |
free_matrix(prop,1,nlstate,iagemin, iagemax+3); |
|
Line 3284 void prevalence(double ***probs, double
|
Line 3506 void prevalence(double ***probs, double
|
| */ |
*/ |
| |
|
| int i, m, jk, j1, bool, z1,j; |
int i, m, jk, j1, bool, z1,j; |
| |
int mi; /* Effective wave */ |
| |
int iage; |
| |
double agebegin, ageend; |
| |
|
| double **prop; |
double **prop; |
| double posprop; |
double posprop; |
|
Line 3303 void prevalence(double ***probs, double
|
Line 3528 void prevalence(double ***probs, double
|
| |
|
| first=1; |
first=1; |
| for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ |
for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ |
| /*for(i1=1; i1<=ncodemax[k1];i1++){ |
for (i=1; i<=nlstate; i++) |
| j1++;*/ |
for(iage=iagemin; iage <= iagemax+3; iage++) |
| |
prop[i][iage]=0.0; |
| for (i=1; i<=nlstate; i++) |
|
| for(m=iagemin; m <= iagemax+3; m++) |
for (i=1; i<=imx; i++) { /* Each individual */ |
| prop[i][m]=0.0; |
bool=1; |
| |
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
| for (i=1; i<=imx; i++) { /* Each individual */ |
for (z1=1; z1<=cptcoveff; z1++) |
| bool=1; |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) |
| if (cptcovn>0) { |
bool=0; |
| for (z1=1; z1<=cptcoveff; z1++) |
} |
| if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) |
if (bool==1) { |
| bool=0; |
/* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ |
| } |
for(mi=1; mi<wav[i];mi++){ |
| if (bool==1) { |
m=mw[mi][i]; |
| for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/ |
agebegin=agev[m][i]; /* Age at beginning of wave before transition*/ |
| |
/* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */ |
| |
if(m >=firstpass && m <=lastpass){ |
| y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ |
y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ |
| if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ |
if ((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 || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); |
if((int)agev[m][i] <iagemin || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); |
| if (s[m][i]>0 && s[m][i]<=nlstate) { |
if (s[m][i]>0 && s[m][i]<=nlstate) { |
| /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ |
/*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ |
| prop[s[m][i]][(int)agev[m][i]] += weight[i]; |
prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ |
| prop[s[m][i]][iagemax+3] += weight[i]; |
prop[s[m][i]][iagemax+3] += weight[i]; |
| } |
} /* end valid statuses */ |
| } |
} /* end selection of dates */ |
| } /* end selection of waves */ |
} /* end selection of waves */ |
| } |
} /* end effective waves */ |
| } |
} /* end bool */ |
| for(i=iagemin; i <= iagemax+3; i++){ |
} |
| for(jk=1,posprop=0; jk <=nlstate ; jk++) { |
for(i=iagemin; i <= iagemax+3; i++){ |
| posprop += prop[jk][i]; |
for(jk=1,posprop=0; jk <=nlstate ; jk++) { |
| } |
posprop += prop[jk][i]; |
| |
} |
| for(jk=1; jk <=nlstate ; jk++){ |
|
| if( i <= iagemax){ |
for(jk=1; jk <=nlstate ; jk++){ |
| if(posprop>=1.e-5){ |
if( i <= iagemax){ |
| probs[i][jk][j1]= prop[jk][i]/posprop; |
if(posprop>=1.e-5){ |
| } else{ |
probs[i][jk][j1]= prop[jk][i]/posprop; |
| if(first==1){ |
} else{ |
| first=0; |
if(first==1){ |
| printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); |
first=0; |
| } |
printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); |
| } |
} |
| } |
} |
| }/* end jk */ |
} |
| }/* end i */ |
}/* end jk */ |
| |
}/* end i */ |
| /*} *//* end i1 */ |
/*} *//* end i1 */ |
| } /* end j1 */ |
} /* end j1 */ |
| |
|
|
Line 3373 void concatwav(int wav[], int **dh, int
|
Line 3601 void concatwav(int wav[], int **dh, int
|
| int i, mi, m; |
int i, mi, m; |
| /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; |
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; |
| double sum=0., jmean=0.;*/ |
double sum=0., jmean=0.;*/ |
| int first; |
int first, firstwo; |
| int j, k=0,jk, ju, jl; |
int j, k=0,jk, ju, jl; |
| double sum=0.; |
double sum=0.; |
| first=0; |
first=0; |
| |
firstwo=0; |
| jmin=100000; |
jmin=100000; |
| jmax=-1; |
jmax=-1; |
| jmean=0.; |
jmean=0.; |
| for(i=1; i<=imx; i++){ |
for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ |
| mi=0; |
mi=0; |
| m=firstpass; |
m=firstpass; |
| while(s[m][i] <= nlstate){ |
while(s[m][i] <= nlstate){ /* a live state */ |
| if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5) |
if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5) |
| mw[++mi][i]=m; |
mw[++mi][i]=m; |
| if(m >=lastpass) |
if(m >=lastpass) |
|
Line 3391 void concatwav(int wav[], int **dh, int
|
Line 3620 void concatwav(int wav[], int **dh, int
|
| else |
else |
| m++; |
m++; |
| }/* end while */ |
}/* end while */ |
| if (s[m][i] > nlstate){ |
if (s[m][i] > nlstate){ /* In a death state */ |
| mi++; /* Death is another wave */ |
mi++; /* Death is another wave */ |
| /* if(mi==0) never been interviewed correctly before death */ |
/* if(mi==0) never been interviewed correctly before death */ |
| /* Only death is a correct wave */ |
/* Only death is a correct wave */ |
| mw[mi][i]=m; |
mw[mi][i]=m; |
| |
}else if (andc[i] != 9999) { /* A death occured after lastpass */ |
| |
m++; |
| |
mi++; |
| |
s[m][i]=nlstate+1; /* We are setting the status to the last of non live state */ |
| |
mw[mi][i]=m; |
| |
nbwarn++; |
| |
if(firstwo==0){ |
| |
printf("Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\nOthers in log file only\n",num[i],i,lastpass,nlstate+1, m); |
| |
fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m); |
| |
firstwo=1; |
| |
} |
| |
if(firstwo==1){ |
| |
fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m); |
| |
} |
| } |
} |
| |
|
| wav[i]=mi; |
wav[i]=mi; |
| if(mi==0){ |
if(mi==0){ |
| nbwarn++; |
nbwarn++; |
|
Line 3410 void concatwav(int wav[], int **dh, int
|
Line 3652 void concatwav(int wav[], int **dh, int
|
| } |
} |
| } /* end mi==0 */ |
} /* end mi==0 */ |
| } /* End individuals */ |
} /* End individuals */ |
| |
/* wav and mw are no more changed */ |
| |
|
| |
|
| for(i=1; i<=imx; i++){ |
for(i=1; i<=imx; i++){ |
| for(mi=1; mi<wav[i];mi++){ |
for(mi=1; mi<wav[i];mi++){ |
| if (stepm <=0) |
if (stepm <=0) |
|
Line 3952 void cvevsij(double ***eij, double x[],
|
Line 4196 void cvevsij(double ***eij, double x[],
|
| } |
} |
| |
|
| /************ Variance ******************/ |
/************ Variance ******************/ |
| void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) |
void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) |
| { |
{ |
| /* Variance of health expectancies */ |
/* Variance of health expectancies */ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ |
|
Line 4006 void cvevsij(double ***eij, double x[],
|
Line 4250 void cvevsij(double ***eij, double x[],
|
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); |
| } |
} |
| printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
| |
|
| fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
| pstamp(ficresprobmorprev); |
pstamp(ficresprobmorprev); |
| fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); |
fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); |
|
Line 4017 void cvevsij(double ***eij, double x[],
|
Line 4260 void cvevsij(double ***eij, double x[],
|
| fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j); |
fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j); |
| } |
} |
| fprintf(ficresprobmorprev,"\n"); |
fprintf(ficresprobmorprev,"\n"); |
| |
|
| fprintf(ficgp,"\n# Routine varevsij"); |
fprintf(ficgp,"\n# Routine varevsij"); |
| fprintf(ficgp,"\nunset title \n"); |
fprintf(ficgp,"\nunset title \n"); |
| /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ |
/* fprintf(fichtm, "#Local time at start: %s", strstart);*/ |
|
Line 4054 void cvevsij(double ***eij, double x[],
|
Line 4298 void cvevsij(double ***eij, double x[],
|
| /* For example we decided to compute the life expectancy with the smallest unit */ |
/* For example we decided to compute the life expectancy with the smallest unit */ |
| /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. |
/* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. |
| nhstepm is the number of hstepm from age to agelim |
nhstepm is the number of hstepm from age to agelim |
| nstepm is the number of stepm from age to agelin. |
nstepm is the number of stepm from age to agelim. |
| Look at function hpijx to understand why (it is linked to memory size questions) */ |
Look at function hpijx to understand why because of memory size limitations, |
| /* We decided (b) to get a life expectancy respecting the most precise curvature of the |
we decided (b) to get a life expectancy respecting the most precise curvature of the |
| survival function given by stepm (the optimization length). Unfortunately it |
survival function given by stepm (the optimization length). Unfortunately it |
| means that if the survival funtion is printed every two years of age and if |
means that if the survival funtion is printed every two years of age and if |
| you sum them up and add 1 year (area under the trapezoids) you won't get the same |
you sum them up and add 1 year (area under the trapezoids) you won't get the same |
|
Line 4077 void cvevsij(double ***eij, double x[],
|
Line 4321 void cvevsij(double ***eij, double x[],
|
| 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); |
| } |
} |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
|
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 4090 void cvevsij(double ***eij, double x[],
|
Line 4334 void cvevsij(double ***eij, double x[],
|
| } |
} |
| } |
} |
| |
|
| |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */ |
| for(j=1; j<= nlstate; j++){ |
for(j=1; j<= nlstate; j++){ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
| gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
| } |
} |
| } |
} |
| /* This for computing probability of death (h=1 means |
/* Next for computing probability of death (h=1 means |
| computed over hstepm matrices product = hstepm*stepm months) |
computed over hstepm matrices product = hstepm*stepm months) |
| as a weighted average of prlim. |
as a weighted average of prlim. |
| */ |
*/ |
|
Line 4108 void cvevsij(double ***eij, double x[],
|
Line 4353 void cvevsij(double ***eij, double x[],
|
| |
|
| 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); |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
|
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear, ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 4121 void cvevsij(double ***eij, double x[],
|
Line 4366 void cvevsij(double ***eij, double x[],
|
| } |
} |
| } |
} |
| |
|
| |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
| |
|
| for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ |
for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
|
Line 4183 void cvevsij(double ***eij, double x[],
|
Line 4430 void cvevsij(double ***eij, double x[],
|
| varppt[j][i]=doldmp[j][i]; |
varppt[j][i]=doldmp[j][i]; |
| /* end ppptj */ |
/* end ppptj */ |
| /* x centered again */ |
/* x centered again */ |
| hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
|
| prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyear,ij); |
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 4200 void cvevsij(double ***eij, double x[],
|
Line 4447 void cvevsij(double ***eij, double x[],
|
| computed over hstepm (estepm) matrices product = hstepm*stepm months) |
computed over hstepm (estepm) matrices product = hstepm*stepm months) |
| as a weighted average of prlim. |
as a weighted average of prlim. |
| */ |
*/ |
| |
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
| for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
| for(i=1,gmp[j]=0.;i<= nlstate; i++) |
for(i=1,gmp[j]=0.;i<= nlstate; i++) |
| gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
|
Line 4262 void cvevsij(double ***eij, double x[],
|
Line 4510 void cvevsij(double ***eij, double x[],
|
| } /* end varevsij */ |
} /* end varevsij */ |
| |
|
| /************ Variance of prevlim ******************/ |
/************ Variance of prevlim ******************/ |
| void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[]) |
void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[]) |
| { |
{ |
| /* Variance of prevalence limit */ |
/* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
| |
|
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
|
Line 4272 void cvevsij(double ***eij, double x[],
|
Line 4520 void cvevsij(double ***eij, double x[],
|
| double *xp; |
double *xp; |
| double *gp, *gm; |
double *gp, *gm; |
| double **gradg, **trgradg; |
double **gradg, **trgradg; |
| |
double **mgm, **mgp; |
| double age,agelim; |
double age,agelim; |
| int theta; |
int theta; |
| |
|
|
Line 4294 void cvevsij(double ***eij, double x[],
|
Line 4543 void cvevsij(double ***eij, double x[],
|
| if (stepm >= YEARM) hstepm=1; |
if (stepm >= YEARM) hstepm=1; |
| nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
| gradg=matrix(1,npar,1,nlstate); |
gradg=matrix(1,npar,1,nlstate); |
| |
mgp=matrix(1,npar,1,nlstate); |
| |
mgm=matrix(1,npar,1,nlstate); |
| gp=vector(1,nlstate); |
gp=vector(1,nlstate); |
| gm=vector(1,nlstate); |
gm=vector(1,nlstate); |
| |
|
|
Line 4301 void cvevsij(double ***eij, double x[],
|
Line 4552 void cvevsij(double ***eij, double x[],
|
| for(i=1; i<=npar; i++){ /* Computes gradient */ |
for(i=1; i<=npar; i++){ /* Computes gradient */ |
| 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,ncvyear,ij); |
if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) |
| for(i=1;i<=nlstate;i++) |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
else |
| |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
for(i=1;i<=nlstate;i++){ |
| gp[i] = prlim[i][i]; |
gp[i] = prlim[i][i]; |
| |
mgp[theta][i] = prlim[i][i]; |
| |
} |
| for(i=1; i<=npar; i++) /* Computes gradient */ |
for(i=1; i<=npar; i++) /* Computes gradient */ |
| 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,ncvyear,ij); |
if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) |
| for(i=1;i<=nlstate;i++) |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
else |
| |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
for(i=1;i<=nlstate;i++){ |
| gm[i] = prlim[i][i]; |
gm[i] = prlim[i][i]; |
| |
mgm[theta][i] = prlim[i][i]; |
| |
} |
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; |
gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; |
| |
/* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */ |
| } /* End theta */ |
} /* End theta */ |
| |
|
| trgradg =matrix(1,nlstate,1,npar); |
trgradg =matrix(1,nlstate,1,npar); |
|
Line 4320 void cvevsij(double ***eij, double x[],
|
Line 4580 void cvevsij(double ***eij, double x[],
|
| for(j=1; j<=nlstate;j++) |
for(j=1; j<=nlstate;j++) |
| for(theta=1; theta <=npar; theta++) |
for(theta=1; theta <=npar; theta++) |
| trgradg[j][theta]=gradg[theta][j]; |
trgradg[j][theta]=gradg[theta][j]; |
| |
/* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ |
| |
/* printf("\nmgm mgp %d ",(int)age); */ |
| |
/* for(j=1; j<=nlstate;j++){ */ |
| |
/* printf(" %d ",j); */ |
| |
/* for(theta=1; theta <=npar; theta++) */ |
| |
/* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */ |
| |
/* printf("\n "); */ |
| |
/* } */ |
| |
/* } */ |
| |
/* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ |
| |
/* printf("\n gradg %d ",(int)age); */ |
| |
/* for(j=1; j<=nlstate;j++){ */ |
| |
/* printf("%d ",j); */ |
| |
/* for(theta=1; theta <=npar; theta++) */ |
| |
/* printf("%d %lf ",theta,gradg[theta][j]); */ |
| |
/* printf("\n "); */ |
| |
/* } */ |
| |
/* } */ |
| |
|
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| varpl[i][(int)age] =0.; |
varpl[i][(int)age] =0.; |
| |
if((int)age==79 ||(int)age== 80 ||(int)age== 81){ |
| |
matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); |
| |
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); |
| |
}else{ |
| matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); |
matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); |
| matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); |
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); |
| |
} |
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ |
varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ |
| |
|
|
Line 4334 void cvevsij(double ***eij, double x[],
|
Line 4617 void cvevsij(double ***eij, double x[],
|
| fprintf(ficresvpl,"\n"); |
fprintf(ficresvpl,"\n"); |
| free_vector(gp,1,nlstate); |
free_vector(gp,1,nlstate); |
| free_vector(gm,1,nlstate); |
free_vector(gm,1,nlstate); |
| |
free_matrix(mgm,1,npar,1,nlstate); |
| |
free_matrix(mgp,1,npar,1,nlstate); |
| free_matrix(gradg,1,npar,1,nlstate); |
free_matrix(gradg,1,npar,1,nlstate); |
| free_matrix(trgradg,1,nlstate,1,npar); |
free_matrix(trgradg,1,nlstate,1,npar); |
| } /* End age */ |
} /* End age */ |
|
Line 4684 To be simple, these graphs help to under
|
Line 4969 To be simple, these graphs help to under
|
| void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ |
void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ |
| int lastpass, int stepm, int weightopt, char model[],\ |
int lastpass, int stepm, int weightopt, char model[],\ |
| int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ |
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ |
| int popforecast, int estepm ,\ |
int popforecast, int prevfcast, int estepm , \ |
| double jprev1, double mprev1,double anprev1, \ |
double jprev1, double mprev1,double anprev1, double dateprev1, \ |
| double jprev2, double mprev2,double anprev2){ |
double jprev2, double mprev2,double anprev2, double dateprev2){ |
| int jj1, k1, i1, cpt; |
int jj1, k1, i1, cpt; |
| |
|
| fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
| <li><a href='#secondorder'>Result files (second order (variance)</a>\n \ |
<li><a href='#secondorder'>Result files (second order (variance)</a>\n \ |
| </ul>"); |
</ul>"); |
| fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \ |
fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n"); |
| - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ", |
fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n", |
| jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_")); |
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm")); |
| |
fprintf(fichtm,"<li> - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file) ", |
| |
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm")); |
| |
fprintf(fichtm,", <a href=\"%s\">%s</a> (text file) <br>\n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", |
- Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", |
| stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_")); |
stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_")); |
|
Line 4702 void printinghtml(char fileresu[], char
|
Line 4990 void printinghtml(char fileresu[], char
|
| - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); |
subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ |
- (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ |
| <a href=\"%s\">%s</a> <br>\n", |
<a href=\"%s\">%s</a> <br>\n", |
| estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); |
estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); |
| fprintf(fichtm,"\ |
if(prevfcast==1){ |
| - Population projections by age and states: \ |
fprintf(fichtm,"\ |
| |
- Prevalence projections by age and states: \ |
| <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); |
<a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); |
| |
} |
| |
|
| fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
| |
|
|
Line 4727 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 5017 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| } |
} |
| /* aij, bij */ |
/* aij, bij */ |
| fprintf(fichtm,"<br>- Logit model, 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.svg\">%s_%d-1.svg</a><br> \ |
| <img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Pij */ |
/* Pij */ |
| fprintf(fichtm,"<br>\n- Pij 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.svg\">%s_%d-2.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.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Quasi-incidences */ |
/* Quasi-incidences */ |
| fprintf(fichtm,"<br>\n- Iij 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 hPij \ |
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: hPij/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.svg\">%s_%d-3.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.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Survival functions (period) in state j */ |
/* Survival functions (period) in state j */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
|
Line 4745 divided by h: hPij/h : <a href=\"%s_%d-3
|
Line 5035 divided by h: hPij/h : <a href=\"%s_%d-3
|
| } |
} |
| /* 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 any different live states 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.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); |
| } |
} |
|
Line 4754 divided by h: hPij/h : <a href=\"%s_%d-3
|
Line 5044 divided by h: hPij/h : <a href=\"%s_%d-3
|
| 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.svg\">%s%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.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); |
| } |
} |
| |
if(prevfcast==1){ |
| |
/* Projection of prevalence up to period (stable) prevalence in each health state */ |
| |
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> \ |
| |
<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); |
| |
} |
| |
} |
| |
|
| 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.svg\">%s_%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.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); |
| } |
} |
| /* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
|
Line 4825 See page 'Matrix of variance-covariance
|
Line 5123 See page 'Matrix of variance-covariance
|
| } |
} |
| for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \ |
fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \ |
| prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg <br>\ |
prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d%d.svg\"> %s_%d-%d.svg <br>\ |
| <img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); |
<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); |
| } |
} |
| 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: %s_%d.svg<br>\ |
observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg<br>\ |
| <img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); |
<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); |
| /* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
| }/* End k1 */ |
}/* End k1 */ |
| fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
|
Line 4841 true period expectancies (those weighted
|
Line 5139 true period expectancies (those weighted
|
| } |
} |
| |
|
| /******************* Gnuplot file **************/ |
/******************* Gnuplot file **************/ |
| void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, char pathc[], double p[]){ |
| |
|
| char dirfileres[132],optfileres[132]; |
char dirfileres[132],optfileres[132]; |
| int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; |
int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; |
| |
int lv=0, vlv=0, kl=0; |
| int ng=0; |
int ng=0; |
| int vpopbased; |
int vpopbased; |
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
|
Line 4862 void printinggnuplot(char fileresu[], ch
|
Line 5161 void printinggnuplot(char fileresu[], ch
|
| fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); |
fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); |
| fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";"); |
fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";"); |
| /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */ |
/* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */ |
| fprintf(ficgp,"\nset ter png size 640, 480"); |
fprintf(ficgp,"\nset ter pngcairo size 640, 480"); |
| /* good for mle=4 plot by number of matrix products. |
/* nice for mle=4 plot by number of matrix products. |
| replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */ |
replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */ |
| /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */ |
/* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */ |
| /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */ |
/* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */ |
| fprintf(ficgp,"\nset out \"%s.png\";",subdirf2(optionfilefiname,"ILK_")); |
fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_")); |
| fprintf(ficgp,"\nplot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); |
fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); |
| fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); |
fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_")); |
| |
fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); |
| |
for (i=1; i<= nlstate ; i ++) { |
| |
fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i); |
| |
fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk)); |
| |
fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1); |
| |
for (j=2; j<= nlstate+ndeath ; j ++) { |
| |
fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j); |
| |
} |
| |
fprintf(ficgp,";\nset out; unset ylabel;\n"); |
| |
} |
| |
/* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */ |
| |
/* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */ |
| |
/* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */ |
| fprintf(ficgp,"\nset out;unset log\n"); |
fprintf(ficgp,"\nset out;unset log\n"); |
| /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */ |
/* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */ |
| |
|
| strcpy(dirfileres,optionfilefiname); |
strcpy(dirfileres,optionfilefiname); |
| strcpy(optfileres,"vpl"); |
strcpy(optfileres,"vpl"); |
| /* 1eme*/ |
/* 1eme*/ |
| fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n"); |
for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ |
| for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */ |
| for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
| |
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); |
| |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1); |
| fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); |
fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \n\ |
|
Line 4905 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5228 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| } /* k1 */ |
} /* k1 */ |
| } /* cpt */ |
} /* cpt */ |
| /*2 eme*/ |
/*2 eme*/ |
| fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); |
|
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| |
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); |
| |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); |
fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); |
| 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) |
|
Line 4939 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5272 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| } /* 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.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ |
| } /* k1 */ |
} /* k1 */ |
| |
|
| |
|
| /*3eme*/ |
/*3eme*/ |
| |
|
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
| |
fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%d state=%d",k1, cpt); |
| |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| /* 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.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); |
|
Line 4968 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5313 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| /* Survival functions (period) from state i in state j by initial state i */ |
/* Survival functions (period) from state i in state j by initial state i */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ |
for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| k=3; |
fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt); |
| fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt); |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); |
| 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\n\ |
set ter svg size 640, 480\n\ |
| unset log y\n\ |
unset log y\n\ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| |
k=3; |
| for (i=1; i<= nlstate ; i ++){ |
for (i=1; i<= nlstate ; i ++){ |
| if(i==1) |
if(i==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
|
Line 4993 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5348 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| /* Survival functions (period) from state i in state j by final state j */ |
/* Survival functions (period) from state i in state j by final state j */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ |
for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ |
| k=3; |
|
| fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt); |
| |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
| set ter svg size 640, 480\n\ |
set ter svg size 640, 480\n\ |
| unset log y\n\ |
unset log y\n\ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| |
k=3; |
| for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ |
for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ |
| if(j==1) |
if(j==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
|
Line 5025 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5390 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| } /* end covariate */ |
} /* end covariate */ |
| |
|
| /* CV preval stable (period) for each covariate */ |
/* CV preval stable (period) for each covariate */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ |
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| k=3; |
fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); |
| fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt); |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
| set ter svg size 640, 480\n\ |
set ter svg size 640, 480\n\ |
| unset log y\n\ |
unset log y\n\ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| |
k=3; /* Offset */ |
| for (i=1; i<= nlstate ; i ++){ |
for (i=1; i<= nlstate ; i ++){ |
| if(i==1) |
if(i==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
|
Line 5049 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5424 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| } /* end cpt state*/ |
} /* end cpt state*/ |
| } /* end covariate */ |
} /* end covariate */ |
| |
|
| |
if(prevfcast==1){ |
| |
/* Projection from cross-sectional to stable (period) for each covariate */ |
| |
|
| |
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ |
| |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| |
fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); |
| |
for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\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,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
| |
set ter svg size 640, 480\n\ |
| |
unset log y\n\ |
| |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| |
for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ |
| |
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| |
/*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| |
if(i==1){ |
| |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); |
| |
}else{ |
| |
fprintf(ficgp,",\\\n '' "); |
| |
} |
| |
if(cptcoveff ==0){ /* No covariate */ |
| |
fprintf(ficgp," u 2:("); /* Age is in 2 */ |
| |
/*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ |
| |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
| |
if(i==nlstate+1) |
| |
fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \ |
| |
2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| |
else |
| |
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ |
| |
2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| |
}else{ |
| |
fprintf(ficgp,"u 6:(("); /* Age is in 6 */ |
| |
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| |
kl=0; |
| |
for (k=1; k<=cptcoveff; k++){ /* For each covariate */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
kl++; |
| |
/* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ |
| |
/*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ |
| |
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
| |
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ |
| |
if(k==cptcoveff) |
| |
if(i==nlstate+1) |
| |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \ |
| |
6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| |
else |
| |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \ |
| |
6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| |
else{ |
| |
fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]); |
| |
kl++; |
| |
} |
| |
} /* end covariate */ |
| |
} /* end if covariate */ |
| |
} /* nlstate */ |
| |
fprintf(ficgp,"\nset out\n"); |
| |
} /* end cpt state*/ |
| |
} /* end covariate */ |
| |
} /* End if prevfcast */ |
| |
|
| |
|
| /* proba elementaires */ |
/* proba elementaires */ |
| fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); |
fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); |
| for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
|
Line 5240 void prevforecast(char fileres[], double
|
Line 5693 void prevforecast(char fileres[], double
|
| char fileresf[FILENAMELENGTH]; |
char fileresf[FILENAMELENGTH]; |
| |
|
| agelim=AGESUP; |
agelim=AGESUP; |
| |
/* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people |
| |
in each health status at the date of interview (if between dateprev1 and dateprev2). |
| |
We still use firstpass and lastpass as another selection. |
| |
*/ |
| |
/* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */ |
| |
/* firstpass, lastpass, stepm, weightopt, model); */ |
| prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
| |
|
| strcpy(fileresf,"F_"); |
strcpy(fileresf,"F_"); |
|
Line 5248 void prevforecast(char fileres[], double
|
Line 5707 void prevforecast(char fileres[], double
|
| printf("Problem with forecast resultfile: %s\n", fileresf); |
printf("Problem with forecast resultfile: %s\n", fileresf); |
| fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf); |
fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf); |
| } |
} |
| printf("Computing forecasting: result on file '%s' \n", fileresf); |
printf("Computing forecasting: result on file '%s', please wait... \n", fileresf); |
| fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf); |
fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf); |
| |
|
| if (cptcoveff==0) ncodemax[cptcoveff]=1; |
if (cptcoveff==0) ncodemax[cptcoveff]=1; |
| |
|
|
Line 5290 void prevforecast(char fileres[], double
|
Line 5749 void prevforecast(char fileres[], double
|
| for(cptcov=1, k=0;cptcov<=i1;cptcov++){ |
for(cptcov=1, k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ |
| k=k+1; |
k=k+1; |
| fprintf(ficresf,"\n#******"); |
fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); |
| for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
| fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| } |
} |
| fprintf(ficresf,"******\n"); |
fprintf(ficresf," yearproj age"); |
| fprintf(ficresf,"# Covariate valuofcovar yearproj age"); |
|
| for(j=1; j<=nlstate+ndeath;j++){ |
for(j=1; j<=nlstate+ndeath;j++){ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| fprintf(ficresf," p%d%d",i,j); |
fprintf(ficresf," p%d%d",i,j); |
|
Line 5345 void prevforecast(char fileres[], double
|
Line 5803 void prevforecast(char fileres[], double
|
| if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| |
|
| fclose(ficresf); |
fclose(ficresf); |
| |
printf("End of Computing forecasting \n"); |
| |
fprintf(ficlog,"End of Computing forecasting\n"); |
| |
|
| } |
} |
| |
|
| /************** Forecasting *****not tested NB*************/ |
/************** Forecasting *****not tested NB*************/ |
|
Line 6184 int calandcheckages(int imx, int maxwav,
|
Line 6645 int calandcheckages(int imx, int maxwav,
|
| for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { |
| agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); |
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); |
| for(m=firstpass; (m<= lastpass); m++){ |
for(m=firstpass; (m<= lastpass); m++){ |
| if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ |
if(s[m][i] >0 || s[m][i]==-1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ /* What if s[m][i]=-1 */ |
| if (s[m][i] >= nlstate+1) { |
if (s[m][i] >= nlstate+1) { |
| if(agedc[i]>0){ |
if(agedc[i]>0){ |
| if((int)moisdc[i]!=99 && (int)andc[i]!=9999){ |
if((int)moisdc[i]!=99 && (int)andc[i]!=9999){ |
| agev[m][i]=agedc[i]; |
agev[m][i]=agedc[i]; |
| /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ |
/*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ |
| }else { |
}else { |
| if ((int)andc[i]!=9999){ |
if ((int)andc[i]!=9999){ |
| nbwarn++; |
nbwarn++; |
|
Line 6199 int calandcheckages(int imx, int maxwav,
|
Line 6660 int calandcheckages(int imx, int maxwav,
|
| } |
} |
| } |
} |
| } /* agedc > 0 */ |
} /* agedc > 0 */ |
| } |
} /* end if */ |
| else if(s[m][i] !=9){ /* Standard case, age in fractional |
else if(s[m][i] !=9){ /* Standard case, age in fractional |
| years but with the precision of a month */ |
years but with the precision of a month */ |
| agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); |
agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); |
|
Line 6215 int calandcheckages(int imx, int maxwav,
|
Line 6676 int calandcheckages(int imx, int maxwav,
|
| } |
} |
| /*agev[m][i]=anint[m][i]-annais[i];*/ |
/*agev[m][i]=anint[m][i]-annais[i];*/ |
| /* agev[m][i] = age[i]+2*m;*/ |
/* agev[m][i] = age[i]+2*m;*/ |
| } |
} /* en if 9*/ |
| else { /* =9 */ |
else { /* =9 */ |
| |
/* printf("Debug num[%d]=%ld s[%d][%d]=%d\n",i,num[i], m,i, s[m][i]); */ |
| agev[m][i]=1; |
agev[m][i]=1; |
| s[m][i]=-1; |
s[m][i]=-1; |
| } |
} |
| } |
} |
| else /*= 0 Unknown */ |
else if(s[m][i]==0) /*= 0 Unknown */ |
| agev[m][i]=1; |
agev[m][i]=1; |
| } |
else{ |
| |
printf("Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]); |
| |
fprintf(ficlog, "Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]); |
| |
agev[m][i]=0; |
| |
} |
| |
} /* End for lastpass */ |
| } |
} |
| |
|
| for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { |
| for(m=firstpass; (m<=lastpass); m++){ |
for(m=firstpass; (m<=lastpass); m++){ |
| if (s[m][i] > (nlstate+ndeath)) { |
if (s[m][i] > (nlstate+ndeath)) { |
|
Line 6451 void syscompilerinfo(int logged)
|
Line 6918 void syscompilerinfo(int logged)
|
| |
|
| } |
} |
| |
|
| int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyear){ |
int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){ |
| /*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
| int i, j, k, i1 ; |
int i, j, k, i1 ; |
| /* double ftolpl = 1.e-10; */ |
/* double ftolpl = 1.e-10; */ |
|
Line 6507 void syscompilerinfo(int logged)
|
Line 6974 void syscompilerinfo(int logged)
|
| |
|
| for (age=agebase; age<=agelim; age++){ |
for (age=agebase; age<=agelim; age++){ |
| /* for (age=agebase; age<=agebase; age++){ */ |
/* for (age=agebase; age<=agebase; age++){ */ |
| prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); |
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k); |
| fprintf(ficrespl,"%.0f ",age ); |
fprintf(ficrespl,"%.0f ",age ); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
|
Line 6516 void syscompilerinfo(int logged)
|
Line 6983 void syscompilerinfo(int logged)
|
| tot += prlim[i][i]; |
tot += prlim[i][i]; |
| fprintf(ficrespl," %.5f", prlim[i][i]); |
fprintf(ficrespl," %.5f", prlim[i][i]); |
| } |
} |
| fprintf(ficrespl," %.3f %d\n", tot, *ncvyear); |
fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp); |
| } /* Age */ |
} /* Age */ |
| /* was end of cptcod */ |
/* was end of cptcod */ |
| } /* cptcov */ |
} /* cptcov */ |
|
Line 6609 int main(int argc, char *argv[])
|
Line 7076 int main(int argc, char *argv[])
|
| #endif |
#endif |
| int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); |
int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); |
| int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; |
int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; |
| int ncvyearnp=0; |
int ncvyear=0; /* Number of years needed for the period prevalence to converge */ |
| int *ncvyear=&ncvyearnp; /* Number of years needed for the period prevalence to converge */ |
|
| int jj, ll, li, lj, lk; |
int jj, ll, li, lj, lk; |
| int numlinepar=0; /* Current linenumber of parameter file */ |
int numlinepar=0; /* Current linenumber of parameter file */ |
| int num_filled; |
int num_filled; |
|
Line 6718 int main(int argc, char *argv[])
|
Line 7184 int main(int argc, char *argv[])
|
| printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); |
printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); |
| if(argc <=1){ |
if(argc <=1){ |
| printf("\nEnter the parameter file name: "); |
printf("\nEnter the parameter file name: "); |
| fgets(pathr,FILENAMELENGTH,stdin); |
if(!fgets(pathr,FILENAMELENGTH,stdin)){ |
| |
printf("ERROR Empty parameter file name\n"); |
| |
goto end; |
| |
} |
| i=strlen(pathr); |
i=strlen(pathr); |
| if(pathr[i-1]=='\n') |
if(pathr[i-1]=='\n') |
| pathr[i-1]='\0'; |
pathr[i-1]='\0'; |
| i=strlen(pathr); |
i=strlen(pathr); |
| if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */ |
if(i >= 1 && pathr[i-1]==' ') {/* This may happen when dragging on oS/X! */ |
| pathr[i-1]='\0'; |
pathr[i-1]='\0'; |
| for (tok = pathr; tok != NULL; ){ |
} |
| |
i=strlen(pathr); |
| |
if( i==0 ){ |
| |
printf("ERROR Empty parameter file name\n"); |
| |
goto end; |
| |
} |
| |
for (tok = pathr; tok != NULL; ){ |
| printf("Pathr |%s|\n",pathr); |
printf("Pathr |%s|\n",pathr); |
| while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); |
while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); |
| printf("val= |%s| pathr=%s\n",val,pathr); |
printf("val= |%s| pathr=%s\n",val,pathr); |
|
Line 6861 int main(int argc, char *argv[])
|
Line 7336 int main(int argc, char *argv[])
|
| if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ |
if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ |
| &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ |
&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ |
| if (num_filled != 8) { |
if (num_filled != 8) { |
| printf("Not 8\n"); |
printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); |
| |
printf("but line=%s\n",line); |
| } |
} |
| printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt); |
printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt); |
| } |
} |
| /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ |
/* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ |
| ftolpl=6.e-3; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ |
/*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ |
| /* Third parameter line */ |
/* Third parameter line */ |
| while(fgets(line, MAXLINE, ficpar)) { |
while(fgets(line, MAXLINE, ficpar)) { |
| /* If line starts with a # it is a comment */ |
/* If line starts with a # it is a comment */ |
|
Line 7215 Please run with mle=-1 to get a correct
|
Line 7691 Please run with mle=-1 to get a correct
|
| free_vector(annais,1,n); |
free_vector(annais,1,n); |
| /* free_matrix(mint,1,maxwav,1,n); |
/* free_matrix(mint,1,maxwav,1,n); |
| free_matrix(anint,1,maxwav,1,n);*/ |
free_matrix(anint,1,maxwav,1,n);*/ |
| free_vector(moisdc,1,n); |
/* free_vector(moisdc,1,n); */ |
| free_vector(andc,1,n); |
/* free_vector(andc,1,n); */ |
| /* */ |
/* */ |
| |
|
| wav=ivector(1,imx); |
wav=ivector(1,imx); |
| dh=imatrix(1,lastpass-firstpass+1,1,imx); |
/* dh=imatrix(1,lastpass-firstpass+1,1,imx); */ |
| bh=imatrix(1,lastpass-firstpass+1,1,imx); |
/* bh=imatrix(1,lastpass-firstpass+1,1,imx); */ |
| mw=imatrix(1,lastpass-firstpass+1,1,imx); |
/* mw=imatrix(1,lastpass-firstpass+1,1,imx); */ |
| |
dh=imatrix(1,lastpass-firstpass+2,1,imx); /* We are adding a wave if status is unknown at last wave but death occurs after last wave.*/ |
| |
bh=imatrix(1,lastpass-firstpass+2,1,imx); |
| |
mw=imatrix(1,lastpass-firstpass+2,1,imx); |
| |
|
| /* Concatenates waves */ |
/* Concatenates waves */ |
| |
/* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i. |
| |
Death is a valid wave (if date is known). |
| |
mw[mi][i] is the number of (mi=1 to wav[i]) effective wave out of mi of individual i |
| |
dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] |
| |
and mw[mi+1][i]. dh depends on stepm. |
| |
*/ |
| |
|
| 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); |
| /* */ |
/* */ |
| |
|
| |
free_vector(moisdc,1,n); |
| |
free_vector(andc,1,n); |
| |
|
| /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ |
/* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ |
| |
|
| nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
|
Line 7235 Please run with mle=-1 to get a correct
|
Line 7724 Please run with mle=-1 to get a correct
|
| Ndum =ivector(-1,NCOVMAX); |
Ndum =ivector(-1,NCOVMAX); |
| if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ |
if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ |
| tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ |
tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ |
| /* Nbcode gives the value of the lth modality of jth covariate, in |
/* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in |
| V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ |
V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ |
| /* 1 to ncodemax[j] is the maximum value of this jth covariate */ |
/* 1 to ncodemax[j] which is the maximum value of this jth covariate */ |
| |
|
| /* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ |
/* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ |
| /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/ |
/*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/ |
| /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ |
/* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ |
| |
/* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j, |
| |
* codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded |
| |
* (currently 0 or 1) in the data. |
| |
* In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of |
| |
* corresponding modality (h,j). |
| |
*/ |
| |
|
| h=0; |
h=0; |
| |
|
| |
|
|
Line 7251 Please run with mle=-1 to get a correct
|
Line 7747 Please run with mle=-1 to get a correct
|
| m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
| |
|
| /**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 |
/**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 |
| * For k=4 covariates, h goes from 1 to 2**k |
* For k=4 covariates, h goes from 1 to m=2**k |
| * codtabm(h,k)= 1 & (h-1) >> (k-1) ; |
* codtabm(h,k)= (1 & (h-1) >> (k-1)) + 1; |
| |
* #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| * h\k 1 2 3 4 |
* h\k 1 2 3 4 |
| *______________________________ |
*______________________________ |
| * 1 i=1 1 i=1 1 i=1 1 i=1 1 |
* 1 i=1 1 i=1 1 i=1 1 i=1 1 |
|
Line 7272 Please run with mle=-1 to get a correct
|
Line 7769 Please run with mle=-1 to get a correct
|
| * 15 i=8 1 2 2 2 |
* 15 i=8 1 2 2 2 |
| * 16 2 2 2 2 |
* 16 2 2 2 2 |
| */ |
*/ |
| |
/* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */ |
| |
/* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4 |
| |
* and the value of each covariate? |
| |
* V1=1, V2=1, V3=2, V4=1 ? |
| |
* h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok. |
| |
* h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st. |
| |
* In order to get the real value in the data, we use nbcode |
| |
* nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0 |
| |
* We are keeping this crazy system in order to be able (in the future?) |
| |
* to have more than 2 values (0 or 1) for a covariate. |
| |
* #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| |
* h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1 |
| |
* bbbbbbbb |
| |
* 76543210 |
| |
* h-1 00000101 (6-1=5) |
| |
*(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift |
| |
* & |
| |
* 1 00000001 (1) |
| |
* 00000001 = 1 & ((h-1) >> (k-1)) |
| |
* +1= 00000010 =2 |
| |
* |
| |
* h=14, k=3 => h'=h-1=13, k'=k-1=2 |
| |
* h' 1101 =2^3+2^2+0x2^1+2^0 |
| |
* >>k' 11 |
| |
* & 00000001 |
| |
* = 00000001 |
| |
* +1 = 00000010=2 = codtabm(14,3) |
| |
* Reverse h=6 and m=16? |
| |
* cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1. |
| |
* for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff) |
| |
* decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1 |
| |
* decodtabm(h,j,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (j-1)) & 1) +1 : -1) |
| |
* V3=decodtabm(14,3,2**4)=2 |
| |
* h'=13 1101 =2^3+2^2+0x2^1+2^0 |
| |
*(h-1) >> (j-1) 0011 =13 >> 2 |
| |
* &1 000000001 |
| |
* = 000000001 |
| |
* +1= 000000010 =2 |
| |
* 2211 |
| |
* V1=1+1, V2=0+1, V3=1+1, V4=1+1 |
| |
* V3=2 |
| |
*/ |
| |
|
| /* /\* for(h=1; h <=100 ;h++){ *\/ */ |
/* /\* for(h=1; h <=100 ;h++){ *\/ */ |
| /* /\* printf("h=%2d ", h); *\/ */ |
/* /\* printf("h=%2d ", h); *\/ */ |
| /* /\* for(k=1; k <=10; k++){ *\/ */ |
/* /\* for(k=1; k <=10; k++){ *\/ */ |
|
Line 7319 Please run with mle=-1 to get a correct
|
Line 7859 Please run with mle=-1 to get a correct
|
| printf("Problem with file %s",optionfilegnuplot); |
printf("Problem with file %s",optionfilegnuplot); |
| } |
} |
| else{ |
else{ |
| fprintf(ficgp,"\n# %s\n", version); |
fprintf(ficgp,"\n# IMaCh-%s\n", version); |
| fprintf(ficgp,"# %s\n", optionfilegnuplot); |
fprintf(ficgp,"# %s\n", optionfilegnuplot); |
| //fprintf(ficgp,"set missing 'NaNq'\n"); |
//fprintf(ficgp,"set missing 'NaNq'\n"); |
| fprintf(ficgp,"set datafile missing 'NaNq'\n"); |
fprintf(ficgp,"set datafile missing 'NaNq'\n"); |
|
Line 7346 Please run with mle=-1 to get a correct
|
Line 7886 Please run with mle=-1 to get a correct
|
| else{ |
else{ |
| fprintf(fichtmcov,"<html><head>\n<title>IMaCh Cov %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
fprintf(fichtmcov,"<html><head>\n<title>IMaCh Cov %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=%s<br>\n",\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
| optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
| } |
} |
| |
|
| fprintf(fichtm,"<html><head>\n<title>IMaCh %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本å¦è¡“振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \ |
| <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=%s<br>\n\ |
<font size=\"2\">IMaCh-%s <br> %s</font> \ |
| |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
| |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n\ |
| \n\ |
\n\ |
| <hr size=\"2\" color=\"#EC5E5E\">\ |
<hr size=\"2\" color=\"#EC5E5E\">\ |
| <ul><li><h4>Parameter files</h4>\n\ |
<ul><li><h4>Parameter files</h4>\n\ |
|
Line 7378 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 7920 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 prints on file fileres'p'. */ |
and prints on file fileres'p'. */ |
| freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart); |
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ |
| |
firstpass, lastpass, stepm, weightopt, model); |
| |
|
| fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
| fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
|
Line 7619 Please run with mle=-1 to get a correct
|
Line 8162 Please run with mle=-1 to get a correct
|
| free_matrix(ximort,1,NDIM,1,NDIM); |
free_matrix(ximort,1,NDIM,1,NDIM); |
| #endif |
#endif |
| } /* Endof if mle==-3 mortality only */ |
} /* Endof if mle==-3 mortality only */ |
| /* Standard maximisation */ |
/* Standard */ |
| else{ /* For mle !=- 3 */ |
else{ /* For mle !=- 3, could be 0 or 1 or 4 etc. */ |
| globpr=0;/* debug */ |
globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */ |
| /* Computes likelihood for initial parameters */ |
/* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */ |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
| printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
| for (k=1; k<=npar;k++) |
for (k=1; k<=npar;k++) |
| printf(" %d %8.5f",k,p[k]); |
printf(" %d %8.5f",k,p[k]); |
| printf("\n"); |
printf("\n"); |
| globpr=1; /* again, to print the contributions */ |
if(mle>=1){ /* Could be 1 or 2, Real Maximization */ |
| |
/* mlikeli uses func not funcone */ |
| |
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); |
| |
} |
| |
if(mle==0) {/* No optimization, will print the likelihoods for the datafile */ |
| |
globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */ |
| |
/* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */ |
| |
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
| |
} |
| |
globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */ |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
| printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
| for (k=1; k<=npar;k++) |
for (k=1; k<=npar;k++) |
| printf(" %d %8.5f",k,p[k]); |
printf(" %d %8.5f",k,p[k]); |
| printf("\n"); |
printf("\n"); |
| if(mle>=1){ /* Could be 1 or 2, Real Maximisation */ |
|
| mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); |
|
| } |
|
| |
|
| /*--------- results files --------------*/ |
/*--------- results files --------------*/ |
| fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); |
fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); |
|
Line 7794 Please run with mle=-1 to get a correct
|
Line 8343 Please run with mle=-1 to get a correct
|
| |
|
| fflush(ficlog); |
fflush(ficlog); |
| fflush(ficres); |
fflush(ficres); |
| |
while(fgets(line, MAXLINE, ficpar)) { |
| while((c=getc(ficpar))=='#' && c!= EOF){ |
/* If line starts with a # it is a comment */ |
| ungetc(c,ficpar); |
if (line[0] == '#') { |
| fgets(line, MAXLINE, ficpar); |
numlinepar++; |
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| } |
fputs(line,ficlog); |
| ungetc(c,ficpar); |
continue; |
| |
}else |
| |
break; |
| |
} |
| |
|
| |
/* while((c=getc(ficpar))=='#' && c!= EOF){ */ |
| |
/* ungetc(c,ficpar); */ |
| |
/* fgets(line, MAXLINE, ficpar); */ |
| |
/* fputs(line,stdout); */ |
| |
/* fputs(line,ficparo); */ |
| |
/* } */ |
| |
/* ungetc(c,ficpar); */ |
| |
|
| estepm=0; |
estepm=0; |
| fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); |
if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){ |
| |
|
| |
if (num_filled != 6) { |
| |
printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n"); |
| |
printf("but line=%s\n",line); |
| |
goto end; |
| |
} |
| |
printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl); |
| |
} |
| |
/* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ |
| |
/*ftolpl=6.e-4;*/ /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ |
| |
|
| |
/* fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); */ |
| if (estepm==0 || estepm < stepm) estepm=stepm; |
if (estepm==0 || estepm < stepm) estepm=stepm; |
| if (fage <= 2) { |
if (fage <= 2) { |
| bage = ageminpar; |
bage = ageminpar; |
|
Line 7812 Please run with mle=-1 to get a correct
|
Line 8384 Please run with mle=-1 to get a correct
|
| } |
} |
| |
|
| fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); |
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); |
| fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl); |
| fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl); |
| |
|
| /* Other stuffs, more or less useful */ |
/* Other stuffs, more or less useful */ |
| while((c=getc(ficpar))=='#' && c!= EOF){ |
while((c=getc(ficpar))=='#' && c!= EOF){ |
|
Line 7876 Please run with mle=-1 to get a correct
|
Line 8448 Please run with mle=-1 to get a correct
|
| This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
| Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
| }else |
}else |
| printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p); |
| |
|
| printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
| model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ |
model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \ |
| jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); |
jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2); |
| |
|
| /*------------ free_vector -------------*/ |
/*------------ free_vector -------------*/ |
| /* chdir(path); */ |
/* chdir(path); */ |
| |
|
| free_ivector(wav,1,imx); |
/* free_ivector(wav,1,imx); */ /* Moved after last prevalence call */ |
| free_imatrix(dh,1,lastpass-firstpass+1,1,imx); |
/* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ |
| free_imatrix(bh,1,lastpass-firstpass+1,1,imx); |
/* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ |
| free_imatrix(mw,1,lastpass-firstpass+1,1,imx); |
/* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ |
| free_lvector(num,1,n); |
free_lvector(num,1,n); |
| free_vector(agedc,1,n); |
free_vector(agedc,1,n); |
| /*free_matrix(covar,0,NCOVMAX,1,n);*/ |
/*free_matrix(covar,0,NCOVMAX,1,n);*/ |
|
Line 7903 Please run with mle=-1 to get a correct
|
Line 8475 Please run with mle=-1 to get a correct
|
| /*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
| /*#include "prevlim.h"*/ /* Use ficrespl, ficlog */ |
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */ |
| prlim=matrix(1,nlstate,1,nlstate); |
prlim=matrix(1,nlstate,1,nlstate); |
| prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, ncvyear); |
prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, &ncvyear); |
| fclose(ficrespl); |
fclose(ficrespl); |
| |
|
| #ifdef FREEEXIT2 |
#ifdef FREEEXIT2 |
|
Line 7948 Please run with mle=-1 to get a correct
|
Line 8520 Please run with mle=-1 to get a correct
|
| /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ |
/* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ |
| ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); |
ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); |
| */ |
*/ |
| |
free_ivector(wav,1,imx); |
| |
free_imatrix(dh,1,lastpass-firstpass+2,1,imx); |
| |
free_imatrix(bh,1,lastpass-firstpass+2,1,imx); |
| |
free_imatrix(mw,1,lastpass-firstpass+2,1,imx); |
| |
|
| |
|
| if (mobilav!=0) { |
if (mobilav!=0) { |
| mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
|
Line 7966 Please run with mle=-1 to get a correct
|
Line 8543 Please run with mle=-1 to get a correct
|
| printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); |
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); |
| fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); |
fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); |
| } |
} |
| printf("Computing Health Expectancies: result on file '%s' \n", filerese); |
printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout); |
| fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); |
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog); |
| /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
| |
|
|
Line 7986 Please run with mle=-1 to get a correct
|
Line 8563 Please run with mle=-1 to get a correct
|
| /*}*/ |
/*}*/ |
| } |
} |
| fclose(ficreseij); |
fclose(ficreseij); |
| |
printf("done evsij\n");fflush(stdout); |
| |
fprintf(ficlog,"done evsij\n");fflush(ficlog); |
| |
|
| /*---------- Health expectancies and variances ------------*/ |
/*---------- Health expectancies and variances ------------*/ |
| |
|
|
Line 7997 Please run with mle=-1 to get a correct
|
Line 8575 Please run with mle=-1 to get a correct
|
| printf("Problem with total LE resultfile: %s\n", filerest);goto end; |
printf("Problem with total LE resultfile: %s\n", filerest);goto end; |
| fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; |
fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; |
| } |
} |
| printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); |
printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout); |
| fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); |
fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog); |
| |
|
| |
|
| strcpy(fileresstde,"STDE_"); |
strcpy(fileresstde,"STDE_"); |
|
Line 8007 Please run with mle=-1 to get a correct
|
Line 8585 Please run with mle=-1 to get a correct
|
| printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); |
printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); |
| fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); |
fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); |
| } |
} |
| printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); |
printf(" Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); |
| fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); |
fprintf(ficlog," Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); |
| |
|
| strcpy(filerescve,"CVE_"); |
strcpy(filerescve,"CVE_"); |
| strcat(filerescve,fileresu); |
strcat(filerescve,fileresu); |
|
Line 8016 Please run with mle=-1 to get a correct
|
Line 8594 Please run with mle=-1 to get a correct
|
| printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); |
printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); |
| fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); |
fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); |
| } |
} |
| printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); |
printf(" Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); |
| fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); |
fprintf(ficlog," Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); |
| |
|
| strcpy(fileresv,"V_"); |
strcpy(fileresv,"V_"); |
| strcat(fileresv,fileresu); |
strcat(fileresv,fileresu); |
|
Line 8025 Please run with mle=-1 to get a correct
|
Line 8603 Please run with mle=-1 to get a correct
|
| printf("Problem with variance resultfile: %s\n", fileresv);exit(0); |
printf("Problem with variance resultfile: %s\n", fileresv);exit(0); |
| fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); |
fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); |
| } |
} |
| printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
printf(" Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(stdout); |
| fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
fprintf(ficlog," Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(ficlog); |
| |
|
| /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
| |
|
| for (k=1; k <= (int) pow(2,cptcoveff); k++){ |
for (k=1; k <= (int) pow(2,cptcoveff); k++){ |
| fprintf(ficrest,"\n#****** "); |
fprintf(ficrest,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;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)]); |
| fprintf(ficrest,"******\n"); |
fprintf(ficrest,"******\n"); |
| |
|
| fprintf(ficresstdeij,"\n#****** "); |
fprintf(ficresstdeij,"\n#****** "); |
| fprintf(ficrescveij,"\n#****** "); |
fprintf(ficrescveij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
| fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| } |
} |
| fprintf(ficresstdeij,"******\n"); |
fprintf(ficresstdeij,"******\n"); |
| fprintf(ficrescveij,"******\n"); |
fprintf(ficrescveij,"******\n"); |
| |
|
| fprintf(ficresvij,"\n#****** "); |
fprintf(ficresvij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| fprintf(ficresvij,"******\n"); |
fprintf(ficresvij,"******\n"); |
| |
|
| eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); |
printf(" cvevsij %d, ",k); |
| /* |
fprintf(ficlog, " cvevsij %d, ",k); |
| */ |
cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); |
| /* goto endfree; */ |
printf(" end cvevsij \n "); |
| |
fprintf(ficlog, " end cvevsij \n "); |
| vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
|
| pstamp(ficrest); |
/* |
| |
*/ |
| |
/* goto endfree; */ |
| for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
|
| oldm=oldms;savm=savms; /* ZZ Segmentation fault */ |
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
| cptcod= 0; /* To be deleted */ |
pstamp(ficrest); |
| varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ |
|
| fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); |
|
| if(vpopbased==1) |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); |
oldm=oldms;savm=savms; /* ZZ Segmentation fault */ |
| else |
cptcod= 0; /* To be deleted */ |
| fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); |
printf("varevsij %d \n",vpopbased); |
| fprintf(ficrest,"# Age popbased mobilav e.. (std) "); |
fprintf(ficlog, "varevsij %d \n",vpopbased); |
| for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ |
| fprintf(ficrest,"\n"); |
fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); |
| /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ |
if(vpopbased==1) |
| epj=vector(1,nlstate+1); |
fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); |
| for(age=bage; age <=fage ;age++){ |
else |
| prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */ |
fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); |
| if (vpopbased==1) { |
fprintf(ficrest,"# Age popbased mobilav e.. (std) "); |
| if(mobilav ==0){ |
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
| for(i=1; i<=nlstate;i++) |
fprintf(ficrest,"\n"); |
| prlim[i][i]=probs[(int)age][i][k]; |
/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ |
| }else{ /* mobilav */ |
epj=vector(1,nlstate+1); |
| for(i=1; i<=nlstate;i++) |
printf("Computing age specific period (stable) prevalences in each health state \n"); |
| prlim[i][i]=mobaverage[(int)age][i][k]; |
fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); |
| } |
for(age=bage; age <=fage ;age++){ |
| } |
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */ |
| |
if (vpopbased==1) { |
| fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); |
if(mobilav ==0){ |
| /* printf(" age %4.0f ",age); */ |
for(i=1; i<=nlstate;i++) |
| for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
prlim[i][i]=probs[(int)age][i][k]; |
| for(i=1, epj[j]=0.;i <=nlstate;i++) { |
}else{ /* mobilav */ |
| epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
for(i=1; i<=nlstate;i++) |
| /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ |
prlim[i][i]=mobaverage[(int)age][i][k]; |
| /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ |
|
| } |
|
| epj[nlstate+1] +=epj[j]; |
|
| } |
} |
| /* printf(" age %4.0f \n",age); */ |
} |
| |
|
| for(i=1, vepp=0.;i <=nlstate;i++) |
fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); |
| for(j=1;j <=nlstate;j++) |
/* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */ |
| vepp += vareij[i][j][(int)age]; |
/* printf(" age %4.0f ",age); */ |
| fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
| for(j=1;j <=nlstate;j++){ |
for(i=1, epj[j]=0.;i <=nlstate;i++) { |
| fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
| |
/*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ |
| |
/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ |
| } |
} |
| fprintf(ficrest,"\n"); |
epj[nlstate+1] +=epj[j]; |
| |
} |
| |
/* printf(" age %4.0f \n",age); */ |
| |
|
| |
for(i=1, vepp=0.;i <=nlstate;i++) |
| |
for(j=1;j <=nlstate;j++) |
| |
vepp += vareij[i][j][(int)age]; |
| |
fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
| |
for(j=1;j <=nlstate;j++){ |
| |
fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
| } |
} |
| |
fprintf(ficrest,"\n"); |
| } |
} |
| free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
} /* End vpopbased */ |
| free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| free_vector(epj,1,nlstate+1); |
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| |
free_vector(epj,1,nlstate+1); |
| |
printf("done \n");fflush(stdout); |
| |
fprintf(ficlog,"done\n");fflush(ficlog); |
| |
|
| /*}*/ |
/*}*/ |
| } |
} /* End k */ |
| free_vector(weight,1,n); |
free_vector(weight,1,n); |
| free_imatrix(Tvard,1,NCOVMAX,1,2); |
free_imatrix(Tvard,1,NCOVMAX,1,2); |
| free_imatrix(s,1,maxwav+1,1,n); |
free_imatrix(s,1,maxwav+1,1,n); |
|
Line 8126 Please run with mle=-1 to get a correct
|
Line 8717 Please run with mle=-1 to get a correct
|
| fclose(ficrescveij); |
fclose(ficrescveij); |
| fclose(ficresvij); |
fclose(ficresvij); |
| fclose(ficrest); |
fclose(ficrest); |
| |
printf("done Health expectancies\n");fflush(stdout); |
| |
fprintf(ficlog,"done Health expectancies\n");fflush(ficlog); |
| fclose(ficpar); |
fclose(ficpar); |
| |
|
| /*------- Variance of period (stable) prevalence------*/ |
/*------- Variance of period (stable) prevalence------*/ |
|
Line 8136 Please run with mle=-1 to get a correct
|
Line 8729 Please run with mle=-1 to get a correct
|
| printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); |
printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); |
| exit(0); |
exit(0); |
| } |
} |
| printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl); |
printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); |
| |
fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); |
| |
|
| /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
|
Line 8149 Please run with mle=-1 to get a correct
|
Line 8743 Please run with mle=-1 to get a correct
|
| |
|
| varpl=matrix(1,nlstate,(int) bage, (int) fage); |
varpl=matrix(1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, strstart); |
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart); |
| free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
| /*}*/ |
/*}*/ |
| } |
} |
| |
|
| fclose(ficresvpl); |
fclose(ficresvpl); |
| |
printf("done variance-covariance of period prevalence\n");fflush(stdout); |
| |
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); |
| |
|
| /*---------- End : free ----------------*/ |
/*---------- End : free ----------------*/ |
| if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
![[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