/* $Id$
$State$
$Log$
+ Revision 1.196 2015/08/18 23:17:52 brouard
+ Summary: 0.98q5
+
Revision 1.195 2015/08/18 16:28:39 brouard
Summary: Adding a hack for testing purpose
#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
-#define codtabm(h,k) 1 & (h-1) >> (k-1) ;
+#define codtabm(h,k) (1 & (h-1) >> (k-1))+1
#define MAXN 20000
#define YEARM 12. /**< Number of months per year */
#define AGESUP 130
/* $State$ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-char copyright[]="August 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[]="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 fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int m,nb;
long *num;
-int firstpass=0, lastpass=4,*cod, *Tage,*cens;
+int firstpass=0, lastpass=4,*cod, *cens;
int *ncodemax; /* ncodemax[j]= Number of modalities of the j th
covariate for which somebody answered excluding
undefined. Usually 2: 0 and 1. */
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+int *Tage;
int *Ndum; /** Freq of modality (tricode */
int **codtab; /**< codtab=imatrix(1,100,1,10); */
int **Tvard, *Tprod, cptcovprod, *Tvaraff;
cov[3]= agefin*agefin;;
for (k=1; k<=cptcovn;k++) {
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
- /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
+ /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
}
/*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
nbcode[Tvar[j]][1]=0;
nbcode[Tvar[j]][2]=1;
nbcode[Tvar[j]][3]=2;
+ To be continued (not working yet).
*/
- ij=0; /* ij is similar to i but can jumps over null modalities */
- for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 to 1*/
- if (Ndum[i] == 0) { /* If at least one individual responded to this modality k */
+ ij=0; /* ij is similar to i but can jump over null modalities */
+ for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
+ if (Ndum[i] == 0) { /* If nobody responded to this modality k */
break;
}
ij++;
- nbcode[Tvar[j]][ij]=i; /* stores the original modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
+ nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
cptcode = ij; /* New max modality for covar j */
} /* end of loop on modality i=-1 to 1 or more */
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
fprintf(fichtm,"\n");
- fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov);
- fprintf(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n\
- file %s<br>\n",optionfilehtmcov);
- fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\
+ fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4></br>this page is important in order to visualize confidence intervals and especially correlation between disability and recovery</li>\n",optionfilehtmcov);
+ fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
+ fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
and drawn. It helps understanding how is the covariance between two incidences.\
They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
fprintf(fichtm,"\
\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \
- - 95%% confidence intervals and T statistics are in the log file.<br>\n", rfileres,rfileres);
+ - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.<br> \
+But because parameters are usually highly correlated (a higher incidence of disability \
+and a higher incidence of recovery can give very close observed transition) it might \
+be very useful to look not only at linear confidence intervals estimated from the \
+variances but at the covariance matrix. And instead of looking at the estimated coefficients \
+(parameters) of the logistic regression, it might be more meaningful to visualize the \
+covariance matrix of the one-step probabilities. \
+See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
ij=1;/* To be checked else nbcode[0][0] wrong */
for(j=3; j <=ncovmodel-nagesqr; j++) {
- if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
- ij++;
+ /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
+ if(ij <=cptcovage) { /* Bug valgrind */
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ ij++;
+ }
}
else
fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
ij=1;
for(j=3; j <=ncovmodel-nagesqr; j++){
- if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
- fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
- ij++;
+ if(ij <=cptcovage) { /* Bug valgrind */
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ ij++;
+ }
}
else
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
int jj, ll, li, lj, lk;
int numlinepar=0; /* Current linenumber of parameter file */
+ int num_filled;
int itimes;
int NDIM=2;
int vpopbased=0;
double ***mobaverage;
char line[MAXLINE];
- char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
+ char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
+
+ char model[MAXLINE], modeltemp[MAXLINE];
char pathr[MAXLINE], pathimach[MAXLINE];
char *tok, *val; /* pathtot */
int firstobs=1, lastobs=10;
goto end;
}
fprintf(ficlog,"Log filename:%s\n",filelog);
- fprintf(ficlog,"\n%s\n%s",version,fullversion);
+ fprintf(ficlog,"Version %s %s",version,fullversion);
fprintf(ficlog,"\nEnter the parameter file name: \n");
fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
path=%s \n\
optionfilext=%s\n\
optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
- syscompilerinfo(0);
+ syscompilerinfo(1);
printf("Local time (at start):%s",strstart);
fprintf(ficlog,"Local time (at start): %s",strstart);
/* Reads comments: lines beginning with '#' */
numlinepar=0;
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
+
+ /* First parameter line */
+ while(fgets(line, MAXLINE, ficpar)) {
+ /* If line starts with a # it is a comment */
+ if (line[0] == '#') {
+ numlinepar++;
+ fputs(line,stdout);
+ fputs(line,ficparo);
+ fputs(line,ficlog);
+ continue;
+ }else
+ break;
+ }
+ if((num_filled=sscanf(line,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", \
+ title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){
+ if (num_filled != 5) {
+ printf("Should be 5 parameters\n");
+ }
numlinepar++;
- fputs(line,stdout);
- fputs(line,ficparo);
- fputs(line,ficlog);
+ printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
+ }
+ /* Second parameter line */
+ while(fgets(line, MAXLINE, ficpar)) {
+ /* If line starts with a # it is a comment */
+ if (line[0] == '#') {
+ numlinepar++;
+ fputs(line,stdout);
+ fputs(line,ficparo);
+ fputs(line,ficlog);
+ continue;
+ }else
+ break;
+ }
+ 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){
+ if (num_filled != 8) {
+ printf("Not 8\n");
+ }
+ 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);
}
- ungetc(c,ficpar);
- fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
- numlinepar=numlinepar+3; /* In general */
- printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+ /* Third parameter line */
+ while(fgets(line, MAXLINE, ficpar)) {
+ /* If line starts with a # it is a comment */
+ if (line[0] == '#') {
+ numlinepar++;
+ fputs(line,stdout);
+ fputs(line,ficparo);
+ fputs(line,ficlog);
+ continue;
+ }else
+ break;
+ }
+ if((num_filled=sscanf(line,"model=1+age%[^.\n]\n", model)) !=EOF){
+ if (num_filled != 1) {
+ printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
+ fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
+ model[0]='\0';
+ goto end;
+ }
+ else{
+ if (model[0]=='+'){
+ for(i=1; i<=strlen(model);i++)
+ modeltemp[i-1]=model[i];
+ }
+ strcpy(model,modeltemp);
+ }
+ printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout);
+ }
+ /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
+ /* numlinepar=numlinepar+3; /\* In general *\/ */
+ /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */
model[strlen(model)-1]='\0';
- fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
- fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+ fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+ fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
fflush(ficlog);
/* if(model[0]=='#'|| model[0]== '\0'){ */
if(model[0]=='#'){
m=pow(2,cptcoveff);
- for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
- for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */
- for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/
- for(cpt=1; cpt <=pow(2,k-1); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */
- h++;
- if (h>m)
- h=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
* codtabm(h,k)= 1 & (h-1) >> (k-1) ;
* 6 2 1 2 1
* 7 i=4 1 2 2 1
* 8 2 2 2 1
- * 9 i=5 1 i=3 1 i=2 1 1
- * 10 2 1 1 1
- * 11 i=6 1 2 1 1
- * 12 2 2 1 1
- * 13 i=7 1 i=4 1 2 1
- * 14 2 1 2 1
- * 15 i=8 1 2 2 1
- * 16 2 2 2 1
+ * 9 i=5 1 i=3 1 i=2 1 2
+ * 10 2 1 1 2
+ * 11 i=6 1 2 1 2
+ * 12 2 2 1 2
+ * 13 i=7 1 i=4 1 2 2
+ * 14 2 1 2 2
+ * 15 i=8 1 2 2 2
+ * 16 2 2 2 2
*/
- codtab[h][k]=j;
- /* codtab[12][3]=1; */
- /*codtab[h][Tvar[k]]=j;*/
- printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);
- }
- }
- }
- }
+ for(h=1; h <=100 ;h++){
+ /* printf("h=%2d ", h); */
+ for(k=1; k <=10; k++){
+ /* printf("k=%d %d ",k,codtabm(h,k)); */
+ codtab[h][k]=codtabm(h,k);
+ }
+ /* printf("\n"); */
+ }
+ /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
+ /* for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/ */
+ /* for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */
+ /* for(cpt=1; cpt <=pow(2,k-1); cpt++){ /\* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 *\/ */
+ /* h++; */
+ /* if (h>m) */
+ /* h=1; */
+ /* codtab[h][k]=j; */
+ /* /\* codtab[12][3]=1; *\/ */
+ /* /\*codtab[h][Tvar[k]]=j;*\/ */
+ /* /\* printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); *\/ */
+ /* } */
+ /* } */
+ /* } */
+ /* } */
/* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]);
codtab[1][2]=1;codtab[2][2]=2; */
- /* for(i=1; i <=m ;i++){
- for(k=1; k <=cptcovn; k++){
- printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
- }
- printf("\n");
- }
- scanf("%d",i);*/
+ /* for(i=1; i <=m ;i++){ */
+ /* for(k=1; k <=cptcovn; k++){ */
+ /* printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+ /* scanf("%d",i);*/
free_ivector(Ndum,-1,NCOVMAX);
ftolhess=ftol; /* Usually correct */
hesscov(matcov, p, npar, delti, ftolhess, func);
}
- printf("Parameters and 95%% confidence intervals\n");
- fprintf(ficlog, "Parameters, T and confidence intervals\n");
+ printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
+ fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
for(i=1,jk=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
printf("%d%d ",i,k);
fprintf(ficlog,"%d%d ",i,k);
for(j=1; j <=ncovmodel; j++){
- printf("%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk]));
- fprintf(ficlog,"%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk]));
+ printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+ fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
jk++;
}
printf("\n");
git://henry.ined.fr / .git / commitdiff