--- imach/src/imach.c 2016/09/15 15:01:13 1.251 +++ imach/src/imach.c 2017/05/13 07:26:12 1.266 @@ -1,6 +1,53 @@ -/* $Id: imach.c,v 1.251 2016/09/15 15:01:13 brouard Exp $ +/* $Id: imach.c,v 1.266 2017/05/13 07:26:12 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.266 2017/05/13 07:26:12 brouard + Summary: Version 0.99r13 (improvements and bugs fixed) + + Revision 1.265 2017/04/26 16:22:11 brouard + Summary: imach 0.99r13 Some bugs fixed + + Revision 1.264 2017/04/26 06:01:29 brouard + Summary: Labels in graphs + + Revision 1.263 2017/04/24 15:23:15 brouard + Summary: to save + + Revision 1.262 2017/04/18 16:48:12 brouard + *** empty log message *** + + Revision 1.261 2017/04/05 10:14:09 brouard + Summary: Bug in E_ as well as in T_ fixed nres-1 vs k1-1 + + Revision 1.260 2017/04/04 17:46:59 brouard + Summary: Gnuplot indexations fixed (humm) + + Revision 1.259 2017/04/04 13:01:16 brouard + Summary: Some errors to warnings only if date of death is unknown but status is death we could set to pi3 + + Revision 1.258 2017/04/03 10:17:47 brouard + Summary: Version 0.99r12 + + Some cleanings, conformed with updated documentation. + + Revision 1.257 2017/03/29 16:53:30 brouard + Summary: Temp + + Revision 1.256 2017/03/27 05:50:23 brouard + Summary: Temporary + + Revision 1.255 2017/03/08 16:02:28 brouard + Summary: IMaCh version 0.99r10 bugs in gnuplot fixed + + Revision 1.254 2017/03/08 07:13:00 brouard + Summary: Fixing data parameter line + + Revision 1.253 2016/12/15 11:59:41 brouard + Summary: 0.99 in progress + + Revision 1.252 2016/09/15 21:15:37 brouard + *** empty log message *** + Revision 1.251 2016/09/15 15:01:13 brouard Summary: not working @@ -128,9 +175,7 @@ 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: Start working on projected prevalences Revision 1.209 2015/11/17 22:12:03 brouard Summary: Adding ftolpl parameter Author: N Brouard @@ -955,12 +1000,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.251 2016/09/15 15:01:13 brouard Exp $ */ +/* $Id: imach.c,v 1.266 2017/05/13 07:26:12 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; -char fullversion[]="$Revision: 1.251 $ $Date: 2016/09/15 15:01:13 $"; +char fullversion[]="$Revision: 1.266 $ $Date: 2017/05/13 07:26:12 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1161,6 +1206,7 @@ int *TvarsQind; #define MAXRESULTLINES 10 int nresult=0; +int parameterline=0; /* # of the parameter (type) line */ int TKresult[MAXRESULTLINES]; int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */ int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */ @@ -2393,8 +2439,8 @@ void powell(double p[], double **xi, int flatd++; } if(flatd >0){ - printf("%d flat directions\n",flatd); - fprintf(ficlog,"%d flat directions\n",flatd); + printf("%d flat directions: ",flatd); + fprintf(ficlog,"%d flat directions :",flatd); for (j=1;j<=n;j++) { if(flatdir[j]>0){ printf("%d ",j); @@ -2578,7 +2624,7 @@ Earliest age to start was %d-%d=%d, ncvl /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */ double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij, int nres) { - /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit + /* Computes the prevalence limit in each live state at age x and for covariate combination ij (<=2**cptcoveff) by left multiplying the unit 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 */ @@ -2614,12 +2660,12 @@ Earliest age to start was %d-%d=%d, ncvl max=vector(1,nlstate); meandiff=vector(1,nlstate); - dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms; - oldm=oldms; savm=savms; - - /* Starting with matrix unity */ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ + dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms; + oldm=oldms; savm=savms; + + /* Starting with matrix unity */ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); } @@ -2639,7 +2685,7 @@ Earliest age to start was %d-%d=%d, ncvl for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; - /* printf("bprevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ + /* printf("bprevalim Dummy agefin=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agefin,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ } /* for (k=1; k<=cptcovn;k++) { */ /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ @@ -2693,8 +2739,27 @@ Earliest age to start was %d-%d=%d, ncvl /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, ageminpar, agemaxpar, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij)); /* Bug Valgrind */ + /* if((int)age == 70){ */ + /* printf(" Backward prevalim age=%d agefin=%d \n", (int) age, (int) agefin); */ + /* for(i=1; i<=nlstate+ndeath; i++) { */ + /* printf("%d newm= ",i); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",newm[i][j]); */ + /* } */ + /* printf("oldm * "); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",oldm[i][j]); */ + /* } */ + /* printf(" pmmij "); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",pmmij[i][j]); */ + /* } */ + /* printf("\n"); */ + /* } */ + /* } */ savm=oldm; oldm=newm; + for(j=1; j<=nlstate; j++){ max[j]=0.; min[j]=1.; @@ -2745,7 +2810,7 @@ Oldest age to start was %d-%d=%d, ncvloo double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) { /* According to parameters values stored in x and the covariate's values stored in cov, - computes the probability to be observed in state j being in state i by appying the + computes the probability to be observed in state j (after stepm years) being in state i by appying the model to the ncovmodel covariates (including constant and age). lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc] and, according on how parameters are entered, the position of the coefficient xij(nc) of the @@ -2754,8 +2819,9 @@ double **pmij(double **ps, double *cov, j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation, sums on j different of i to get 1-pii/pii, deduces pii, and then all pij. - Outputs ps[i][j] the probability to be observed in j being in j according to + Outputs ps[i][j] or probability to be observed in j being in i according to the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij] + Sum on j ps[i][j] should equal to 1. */ double s1, lnpijopii; /*double t34;*/ @@ -2819,7 +2885,7 @@ double **pmij(double **ps, double *cov, /* for(i=1; i<= npar; i++) printf("%f ",x[i]); goto end;*/ - return ps; + return ps; /* Pointer is unchanged since its call */ } /*************** backward transition probabilities ***************/ @@ -2828,8 +2894,8 @@ double **pmij(double **ps, double *cov, /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) { - /* Computes the backward probability at age agefin and covariate ij - * and returns in **ps as well as **bmij. + /* Computes the backward probability at age agefin and covariate combination ij. In fact cov is already filled and x too. + * Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij. */ int i, ii, j,k; @@ -2846,43 +2912,56 @@ double **pmij(double **ps, double *cov, agefin=cov[2]; /* bmij *//* age is cov[2], ij is included in cov, but we need for - the observed prevalence (with this covariate ij) */ - dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); + the observed prevalence (with this covariate ij) at beginning of transition */ + /* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ + pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */ + /* outputs pmmij which is a stochastic matrix */ /* We do have the matrix Px in savm and we need pij */ for (j=1;j<=nlstate+ndeath;j++){ - sumnew=0.; /* w1 p11 + w2 p21 only on live states */ + sumnew=0.; /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */ for (ii=1;ii<=nlstate;ii++){ - sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; + /* sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; */ + sumnew+=pmmij[ii][j]*prevacurrent[(int)agefin][ii][ij]; /* Yes prevalence at beginning of transition */ } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ - for (ii=1;ii<=nlstate+ndeath;ii++){ - if(sumnew >= 1.e-10){ + if(sumnew >= 1.e-10){ + for (ii=1;ii<=nlstate+ndeath;ii++){ /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ /* }else if(agefin >= agemaxpar+stepm/YEARM){ */ /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ /* }else */ doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); - }else{ - ; - /* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */ - } - } /*End ii */ - } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ - /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */ - bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ + } /*End ii */ + }else{ /* We put the identity matrix */ + for (ii=1;ii<=nlstate+ndeath;ii++){ + doldm[ii][j]=(ii==j ? 1. : 0.0); + } /*End ii */ + /* printf("Problem internal bmij A: sum_i w_i*p_ij=N.j/N.. <1.e-10 i=%d, j=%d, sumnew=%lf,agefin=%d\n",ii,j,sumnew, (int)agefin); */ + } + } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] or identity*/ + /* left Product of this diag matrix by dsavm=Px (dnewm=dsavm*doldm) */ + /* bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /\* Bug Valgrind *\/ */ + bbmij=matprod2(dnewm, pmmij,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ /* left Product of this matrix by diag matrix of prevalences (savm) */ for (j=1;j<=nlstate+ndeath;j++){ + sumnew=0.; for (ii=1;ii<=nlstate+ndeath;ii++){ + sumnew+=prevacurrent[(int)agefin][ii][ij]; dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); } - } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */ - ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ + /* if(sumnew <0.9){ */ + /* printf("Problem internal bmij B: sum on i wi <0.9: j=%d, sum_i wi=%lf,agefin=%d\n",j,sumnew, (int)agefin); */ + /* } */ + } /* End j, At the end dsavm is diag[(w_i)] */ + /* What if dsavm doesn't sum ii to 1? */ + /* ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /\* Bug Valgrind *\/ */ + ps=matprod2(ps, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ /* end bmij */ - return ps; + return ps; /*pointer is unchanged */ } /*************** transition probabilities ***************/ @@ -3108,7 +3187,7 @@ double ***hpxij(double ***po, int nhstep /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) { - /* Computes the transition matrix starting at age 'age' over + /* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' months (i.e. until age (in years) age+nhstepm*hstepm*stepm/12) by multiplying nhstepm*hstepm matrices. @@ -3116,18 +3195,19 @@ double ***hbxij(double ***po, int nhstep (typically every 2 years instead of every month which is too big for the memory). Model is determined by parameters x and covariates have to be - included manually here. - + included manually here. Then we use a call to bmij(x and cov) + The addresss of po (p3mat allocated to the dimension of nhstepm) should be stored for output */ int i, j, d, h, k; - double **out, cov[NCOVMAX+1]; - double **newm; + double **out, cov[NCOVMAX+1], **bmij(); + double **newm, ***newmm; double agexact; double agebegin, ageend; double **oldm, **savm; - oldm=oldms;savm=savms; + newmm=po; /* To be saved */ + oldm=oldms;savm=savms; /* Global pointers */ /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) for (j=1;j<=nlstate+ndeath;j++){ @@ -3140,14 +3220,18 @@ double ***hbxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ + agexact=age-((h-1)*hstepm + (d))*stepm/YEARM; /* age just before transition, d or d-1? */ /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ cov[2]=agexact; if(nagesqr==1) cov[3]= agexact*agexact; - for (k=1; k<=cptcovn;k++) - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; - /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + for (k=1; k<=cptcovn;k++){ + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ + /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; + /* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ + + } for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; @@ -3156,11 +3240,10 @@ double ***hbxij(double ***po, int nhstep cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ - /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ /* Careful transposed matrix */ - /* age is in cov[2] */ + /* age is in cov[2], prevacurrent at beginning of transition. */ /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ /* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ @@ -4122,7 +4205,16 @@ void ludcmp(double **a, int n, int *indx big=0.0; for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; - if (big == 0.0) nrerror("Singular matrix in routine ludcmp"); + if (big == 0.0){ + printf(" Singular Hessian matrix at row %d:\n",i); + for (j=1;j<=n;j++) { + printf(" a[%d][%d]=%f,",i,j,a[i][j]); + fprintf(ficlog," a[%d][%d]=%f,",i,j,a[i][j]); + } + fflush(ficlog); + fclose(ficlog); + nrerror("Singular matrix in routine ludcmp"); + } vv[i]=1.0/big; } for (j=1;j<=n;j++) { @@ -4188,17 +4280,84 @@ void pstamp(FILE *fichier) fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); } +int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) { + + /* y=a+bx regression */ + double sumx = 0.0; /* sum of x */ + double sumx2 = 0.0; /* sum of x**2 */ + double sumxy = 0.0; /* sum of x * y */ + double sumy = 0.0; /* sum of y */ + double sumy2 = 0.0; /* sum of y**2 */ + double sume2; /* sum of square or residuals */ + double yhat; + + double denom=0; + int i; + int ne=*no; + + for ( i=ifi, ne=0;i<=ila;i++) { + if(!isfinite(x[i]) || !isfinite(y[i])){ + /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */ + continue; + } + ne=ne+1; + sumx += x[i]; + sumx2 += x[i]*x[i]; + sumxy += x[i] * y[i]; + sumy += y[i]; + sumy2 += y[i]*y[i]; + denom = (ne * sumx2 - sumx*sumx); + /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */ + } + + denom = (ne * sumx2 - sumx*sumx); + if (denom == 0) { + // vertical, slope m is infinity + *b = INFINITY; + *a = 0; + if (r) *r = 0; + return 1; + } + + *b = (ne * sumxy - sumx * sumy) / denom; + *a = (sumy * sumx2 - sumx * sumxy) / denom; + if (r!=NULL) { + *r = (sumxy - sumx * sumy / ne) / /* compute correlation coeff */ + sqrt((sumx2 - sumx*sumx/ne) * + (sumy2 - sumy*sumy/ne)); + } + *no=ne; + for ( i=ifi, ne=0;i<=ila;i++) { + if(!isfinite(x[i]) || !isfinite(y[i])){ + /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */ + continue; + } + ne=ne+1; + yhat = y[i] - *a -*b* x[i]; + sume2 += yhat * yhat ; + + denom = (ne * sumx2 - sumx*sumx); + /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */ + } + *sb = sqrt(sume2/(ne-2)/(sumx2 - sumx * sumx /ne)); + *sa= *sb * sqrt(sumx2/ne); + + return 0; +} + /************ Frequencies ********************/ void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \ int firstpass, int lastpass, int stepm, int weightopt, char model[]) { /* Some frequencies as well as proposing some starting values */ - int i, m, jk, j1, bool, z1,j, k, iv, jj=0; + int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0, s1=1, s2=1; int iind=0, iage=0; int mi; /* Effective wave */ int first; double ***freq; /* Frequencies */ + double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */ + int no; double *meanq; double **meanqt; double *pp, **prop, *posprop, *pospropt; @@ -4251,6 +4410,8 @@ Title=%s
Datafile=%s Firstpass=%d La } fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate)

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); + y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); + x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+4+AGEMARGE); j1=0; @@ -4268,16 +4429,49 @@ Title=%s
Datafile=%s Firstpass=%d La dateintsum=0; k2cpt=0; - for (j = 0; j <= cptcoveff; j+=cptcoveff){ /* j= 0 constant model */ + if(cptcoveff == 0 ) + nl=1; /* Constant and age model only */ + else + nl=2; + + /* if a constant only model, one pass to compute frequency tables and to write it on ficresp */ + /* Loop on nj=1 or 2 if dummy covariates j!=0 + * Loop on j1(1 to 2**cptcoveff) covariate combination + * freq[s1][s2][iage] =0. + * Loop on iind + * ++freq[s1][s2][iage] weighted + * end iind + * if covariate and j!0 + * headers Variable on one line + * endif cov j!=0 + * header of frequency table by age + * Loop on age + * pp[s1]+=freq[s1][s2][iage] weighted + * pos+=freq[s1][s2][iage] weighted + * Loop on s1 initial state + * fprintf(ficresp + * end s1 + * end age + * if j!=0 computes starting values + * end compute starting values + * end j1 + * end nl + */ + for (nj = 1; nj <= nl; nj++){ /* nj= 1 constant model, nl number of loops. */ + if(nj==1) + j=0; /* First pass for the constant */ + else{ + j=cptcoveff; /* Other passes for the covariate values */ + } first=1; - for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */ + for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all covariates combination of the model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */ posproptt=0.; /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ for (i=-5; i<=nlstate+ndeath; i++) - for (jk=-5; jk<=nlstate+ndeath; jk++) + for (s2=-5; s2<=nlstate+ndeath; s2++) for(m=iagemin; m <= iagemax+3; m++) - freq[i][jk][m]=0; + freq[i][s2][m]=0; for (i=1; i<=nlstate; i++) { for(m=iagemin; m <= iagemax+3; m++) @@ -4295,7 +4489,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* dateintsum=0; */ /* k2cpt=0; */ - /* For that combination of covariate j1, we count and print the frequencies in one pass */ + /* For that combination of covariates j1 (V4=1 V3=0 for example), we count and print the frequencies in one pass */ for (iind=1; iind<=imx; iind++) { /* For each individual iind */ bool=1; if(j !=0){ @@ -4311,7 +4505,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* /\* sumnew+=coqvar[z1][iind]; *\/ */ /* }else */ if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */ - /* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */ + /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */ bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), @@ -4322,7 +4516,7 @@ Title=%s
Datafile=%s Firstpass=%d La } /* cptcovn > 0 */ } /* end any */ }/* end j==0 */ - if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */ + if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ /* for(m=firstpass; m<=lastpass; m++){ */ for(mi=1; miDatafile=%s Firstpass=%d La /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - pstamp(ficresp); + if(cptcoveff==0 && nj==1) /* no covariate and first pass */ + pstamp(ficresp); if (cptcoveff>0 && j!=0){ + pstamp(ficresp); printf( "\n#********** Variable "); fprintf(ficresp, "\n#********** Variable "); fprintf(ficresphtm, "\n

********** Variable "); @@ -4413,20 +4609,23 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficlog, "**********\n"); } fprintf(ficresphtm,""); + if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ + fprintf(ficresp, " Age"); + if(nj==2) for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); for(i=1; i<=nlstate;i++) { - fprintf(ficresp, " Age Prev(%d) N(%d) N ",i,i); + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i); fprintf(ficresphtm, "",i,i); } - fprintf(ficresp, "\n"); + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); fprintf(ficresphtm, "\n"); /* Header of frequency table by age */ fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); fprintf(ficresphtmfr," "); - for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(s2=-1; s2 <=nlstate+ndeath; s2++){ for(m=-1; m <=nlstate+ndeath; m++){ - if(jk!=0 && m!=0) - fprintf(ficresphtmfr," ",jk,m); + if(s2!=0 && m!=0) + fprintf(ficresphtmfr," ",s2,m); } } fprintf(ficresphtmfr, "\n"); @@ -4451,94 +4650,112 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtmfr," ",iage); fprintf(ficlog,"Age %d", iage); } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][iage]; + for(s1=1; s1 <=nlstate ; s1++){ + for(m=-1, pp[s1]=0; m <=nlstate+ndeath ; m++) + pp[s1] += freq[s1][m][iage]; } - for(jk=1; jk <=nlstate ; jk++){ + for(s1=1; s1 <=nlstate ; s1++){ for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][iage]; - if(pp[jk]>=1.e-10){ + pos += freq[s1][m][iage]; + if(pp[s1]>=1.e-10){ if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + printf(" %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); } - fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); }else{ if(first==1) - printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); + printf(" %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); + fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); } } - for(jk=1; jk <=nlstate ; jk++){ - /* posprop[jk]=0; */ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ - pp[jk] += freq[jk][m][iage]; - } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ - - for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ - pos += pp[jk]; /* pos is the total number of transitions until this age */ - posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ - pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + for(s1=1; s1 <=nlstate ; s1++){ + /* posprop[s1]=0; */ + for(m=0, pp[s1]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ + pp[s1] += freq[s1][m][iage]; + } /* pp[s1] is the total number of transitions starting from state s1 and any ending status until this age */ + + for(s1=1,pos=0, pospropta=0.; s1 <=nlstate ; s1++){ + pos += pp[s1]; /* pos is the total number of transitions until this age */ + posprop[s1] += prop[s1][iage]; /* prop is the number of transitions from a live state + from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pospropta += prop[s1][iage]; /* prop is the number of transitions from a live state + from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + } + + /* Writing ficresp */ + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + if( iage <= iagemax){ + fprintf(ficresp," %d",iage); + } + }else if( nj==2){ + if( iage <= iagemax){ + fprintf(ficresp," %d",iage); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + } } - for(jk=1; jk <=nlstate ; jk++){ + for(s1=1; s1 <=nlstate ; s1++){ if(pos>=1.e-5){ if(first==1) - printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + printf(" %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos); + fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos); }else{ if(first==1) - printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); + printf(" %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); + fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); } if( iage <= iagemax){ if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - /*probs[iage][jk][j1]= pp[jk]/pos;*/ - /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ - } - else{ - fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage, prop[jk][iage],pospropta); + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + }else if( nj==2){ + fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + } + fprintf(ficresphtm,"",iage,prop[s1][iage]/pospropta, prop[s1][iage],pospropta); + /*probs[iage][s1][j1]= pp[s1]/pos;*/ + /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/ + } else{ + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); + fprintf(ficresphtm,"",iage, prop[s1][iage],pospropta); } } - pospropt[jk] +=posprop[jk]; - } /* end loop jk */ + pospropt[s1] +=posprop[s1]; + } /* end loop s1 */ /* pospropt=0.; */ - for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(s1=-1; s1 <=nlstate+ndeath; s1++){ for(m=-1; m <=nlstate+ndeath; m++){ - if(freq[jk][m][iage] !=0 ) { /* minimizing output */ + if(freq[s1][m][iage] !=0 ) { /* minimizing output */ if(first==1){ - printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); + printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); } - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); + /* printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); */ + fprintf(ficlog," %d%d=%.0f",s1,m,freq[s1][m][iage]); } - if(jk!=0 && m!=0) - fprintf(ficresphtmfr," ",freq[jk][m][iage]); + if(s1!=0 && m!=0) + fprintf(ficresphtmfr," ",freq[s1][m][iage]); } - } /* end loop jk */ + } /* end loop s1 */ posproptt=0.; - for(jk=1; jk <=nlstate; jk++){ - posproptt += pospropt[jk]; + for(s1=1; s1 <=nlstate; s1++){ + posproptt += pospropt[s1]; } fprintf(ficresphtmfr,"\n "); - if(iage <= iagemax){ - fprintf(ficresp,"\n"); - fprintf(ficresphtm,"\n"); + fprintf(ficresphtm,"\n"); + if((cptcoveff==0 && nj==1)|| nj==2 ) { + if(iage <= iagemax) + fprintf(ficresp,"\n"); } if(first==1) printf("Others in log...\n"); fprintf(ficlog,"\n"); } /* end loop age iage */ + fprintf(ficresphtm,""); - for(jk=1; jk <=nlstate ; jk++){ + for(s1=1; s1 <=nlstate ; s1++){ if(posproptt < 1.e-5){ - fprintf(ficresphtm,"",pospropt[jk],posproptt); + fprintf(ficresphtm,"",pospropt[s1],posproptt); }else{ - fprintf(ficresphtm,"",pospropt[jk]/posproptt,pospropt[jk],posproptt); + fprintf(ficresphtm,"",pospropt[s1]/posproptt,pospropt[s1],posproptt); } } fprintf(ficresphtm,"\n"); @@ -4547,7 +4764,8 @@ Title=%s
Datafile=%s Firstpass=%d La if(posproptt < 1.e-5){ fprintf(ficresphtm,"\n

This combination (%d) is not valid and no result will be produced

",j1); fprintf(ficresphtmfr,"\n

This combination (%d) is not valid and no result will be produced

",j1); - fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); + fprintf(ficlog,"# This combination (%d) is not valid and no result will be produced\n",j1); + printf("# This combination (%d) is not valid and no result will be produced\n",j1); invalidvarcomb[j1]=1; }else{ fprintf(ficresphtm,"\n

This combination (%d) is valid and result will be produced.

",j1); @@ -4557,55 +4775,66 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficlog,"\n"); if(j!=0){ printf("#Freqsummary: Starting values for combination j1=%d:\n", j1); - for(i=1,jk=1; i <=nlstate; i++){ + for(i=1,s1=1; i <=nlstate; i++){ for(k=1; k <=(nlstate+ndeath); k++){ if (k != i) { - printf("%d%d ",i,k); - fprintf(ficlog,"%d%d ",i,k); - for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */ - if(jj==1){ /* Constant case */ + for(jj=1; jj <=ncovmodel; jj++){ /* For counting s1 */ + if(jj==1){ /* Constant case (in fact cste + age) */ if(j1==1){ /* All dummy covariates to zero */ freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */ freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */ + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3])); + fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); } - printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3])); - fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); - pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); - }else if( (log(j1-1)/log(2)+1 == jj -2 -nagesqr) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ - pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])); - printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4])); - }else if(jj==2 || nagesqr==1){ /* age or age*age parameter */ - ; + }else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */ + for(iage=iagemin; iage <= iagemax+3; iage++){ + x[iage]= (double)iage; + y[iage]= log(freq[i][k][iage]/freq[i][i][iage]); + /* printf("i=%d, k=%d, s1=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,s1,j1,jj, iage, y[iage]); */ + } + linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */ + pstart[s1]=b; + pstart[s1-1]=a; + }else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ + printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); + printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); + pstart[s1]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])); + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + printf("s1=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",s1,i,k,s1,p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4])); }else{ /* Other cases, like quantitative fixed or varying covariates */ ; } /* printf("%12.7f )", param[i][jj][k]); */ /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ - jk++; + s1++; } /* end jj */ - printf("\n"); - fprintf(ficlog,"\n"); } /* end k!= i */ } /* end k */ - } /* end i, jk */ + } /* end i, s1 */ } /* end j !=0 */ } /* end selected combination of covariate j1 */ if(j==0){ /* We can estimate starting values from the occurences in each case */ printf("#Freqsummary: Starting values for the constants:\n"); fprintf(ficlog,"\n"); - for(i=1,jk=1; i <=nlstate; i++){ + for(i=1,s1=1; i <=nlstate; i++){ for(k=1; k <=(nlstate+ndeath); k++){ if (k != i) { printf("%d%d ",i,k); fprintf(ficlog,"%d%d ",i,k); for(jj=1; jj <=ncovmodel; jj++){ - if(jj==1){ - printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); - fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + pstart[s1]=p[s1]; /* Setting pstart to p values by default */ + if(jj==1){ /* Age has to be done */ + pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); + printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); + fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); } /* printf("%12.7f )", param[i][jj][k]); */ /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ - jk++; + s1++; } printf("\n"); fprintf(ficlog,"\n"); @@ -4614,22 +4843,44 @@ Title=%s
Datafile=%s Firstpass=%d La } printf("#Freqsummary\n"); fprintf(ficlog,"\n"); - for(jk=-1; jk <=nlstate+ndeath; jk++){ - for(m=-1; m <=nlstate+ndeath; m++){ - /* param[i]|j][k]= freq[jk][m][iagemax+3] */ - printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); - /* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */ - /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ - /* fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ + for(s1=-1; s1 <=nlstate+ndeath; s1++){ + for(s2=-1; s2 <=nlstate+ndeath; s2++){ + /* param[i]|j][k]= freq[s1][s2][iagemax+3] */ + printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); + fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); + /* if(freq[s1][s2][iage] !=0 ) { /\* minimizing output *\/ */ + /* printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ + /* fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ /* } */ } - } /* end loop jk */ + } /* end loop s1 */ printf("\n"); fprintf(ficlog,"\n"); } /* end j=0 */ } /* end j */ + + if(mle == -2){ /* We want to use these values as starting values */ + for(i=1, jk=1; i <=nlstate; i++){ + for(j=1; j <=nlstate+ndeath; j++){ + if(j!=i){ + /*ca[0]= k+'a'-1;ca[1]='\0';*/ + printf("%1d%1d",i,j); + fprintf(ficparo,"%1d%1d",i,j); + for(k=1; k<=ncovmodel;k++){ + /* printf(" %lf",param[i][j][k]); */ + /* fprintf(ficparo," %lf",param[i][j][k]); */ + p[jk]=pstart[jk]; + printf(" %f ",pstart[jk]); + fprintf(ficparo," %f ",pstart[jk]); + jk++; + } + printf("\n"); + fprintf(ficparo,"\n"); + } + } + } + } /* end mle=-2 */ dateintmean=dateintsum/k2cpt; fclose(ficresp); @@ -4637,6 +4888,8 @@ Title=%s
Datafile=%s Firstpass=%d La fclose(ficresphtmfr); free_vector(meanq,1,nqfveff); free_matrix(meanqt,1,lastpass,1,nqtveff); + free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); + free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_vector(pospropt,1,nlstate); free_vector(posprop,1,nlstate); @@ -4734,7 +4987,10 @@ void prevalence(double ***probs, double } else{ if(first==1){ first=0; - printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,j1,probs[i][jk][j1]); + printf("Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); + fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); + }else{ + fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); } } } @@ -4794,10 +5050,10 @@ void concatwav(int wav[], int **dh, int #else if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ if(firsthree == 0){ - printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); firsthree=1; } - fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); mw[++mi][i]=m; mli=m; } @@ -4824,13 +5080,9 @@ void concatwav(int wav[], int **dh, int /* if(mi==0) never been interviewed correctly before death */ /* Only death is a correct wave */ mw[mi][i]=m; - } + } /* else not in a death state */ #ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE - else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */ - /* m++; */ - /* mi++; */ - /* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */ - /* mw[mi][i]=m; */ + else if ((int) andc[i] != 9999) { /* Date of death is known */ if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */ nbwarn++; @@ -4843,12 +5095,12 @@ void concatwav(int wav[], int **dh, int }else{ /* Death occured afer last wave potential bias */ nberr++; if(firstwo==0){ - printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); firstwo=1; } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); } - }else{ /* end date of interview is known */ + }else{ /* if date of interview is unknown */ /* death is known but not confirmed by death status at any wave */ if(firstfour==0){ printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); @@ -4879,7 +5131,7 @@ void concatwav(int wav[], int **dh, int if (stepm <=0) dh[mi][i]=1; else{ - if (s[mw[mi+1][i]][i] > nlstate) { /* A death */ + if (s[mw[mi+1][i]][i] > nlstate) { /* A death, but what if date is unknown? */ if (agedc[i] < 2*AGESUP) { j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); if(j==0) j=1; /* Survives at least one month after exam */ @@ -5976,7 +6228,7 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.
    • \n",optionfilehtmcov); + fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. %s
    • \n",optionfilehtmcov,optionfilehtmcov); fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ and drawn. It helps understanding how is the covariance between two incidences.\ @@ -6193,7 +6445,7 @@ To be simple, these graphs help to under fprintf(ficgp,"\nset parametric;unset label"); fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); fprintf(ficgp,"\nset ter svg size 640, 480"); - fprintf(fichtmcov,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ + fprintf(fichtmcov,"\n


    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ : \ %s_%d%1d%1d-%1d%1d.svg, ",k1,l1,k2,l2,\ subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \ @@ -6204,16 +6456,16 @@ To be simple, these graphs help to under fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ - mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ - mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); + mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \ + mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */ }else{ first=0; fprintf(fichtmcov," %d (%.3f),",(int) age, c12); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ - mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ - mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); + mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \ + mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); }/* if first */ } /* age mod 5 */ } /* end loop age */ @@ -6241,7 +6493,7 @@ To be simple, these graphs help to under void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ int lastpass, int stepm, int weightopt, char model[],\ int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ - int popforecast, int prevfcast, int backcast, int estepm , \ + int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \ double jprev1, double mprev1,double anprev1, double dateprev1, \ double jprev2, double mprev2,double anprev2, double dateprev2){ int jj1, k1, i1, cpt, k4, nres; @@ -6279,21 +6531,66 @@ void printinghtml(char fileresu[], char %s
    \n", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); } - fprintf(fichtm," \n

    Age%d%d%d%d
    %d%d%.5f%.0f%.0f%dNaNq%.0f%.0f%d%.5f%.0f%.0f%dNaNq%.0f%.0f%.0f%.0f
    TotNanq%.0f%.0fNanq%.0f%.0f%.5f%.0f%.0f%.5f%.0f%.0f