--- imach/src/imach.c 2002/05/21 18:44:41 1.42 +++ imach/src/imach.c 2015/12/16 06:57:54 1.214 @@ -1,3558 +1,8871 @@ -/* $Id: imach.c,v 1.42 2002/05/21 18:44:41 brouard Exp $ - Interpolated Markov Chain - - Short summary of the programme: - - This program computes Healthy Life Expectancies from - cross-longitudinal data. Cross-longitudinal data consist in: -1- a - first survey ("cross") where individuals from different ages are - interviewed on their health status or degree of disability (in the - case of a health survey which is our main interest) -2- at least a - second wave of interviews ("longitudinal") which measure each change - (if any) in individual health status. Health expectancies are - computed from the time spent in each health state according to a - model. More health states you consider, more time is necessary to reach the - Maximum Likelihood of the parameters involved in the model. The - simplest model is the multinomial logistic model where pij is the - probability to be observed in state j at the second wave - conditional to be observed in state i at the first wave. Therefore - the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where - 'age' is age and 'sex' is a covariate. If you want to have a more - complex model than "constant and age", you should modify the program - where the markup *Covariates have to be included here again* invites - you to do it. More covariates you add, slower the - convergence. - - The advantage of this computer programme, compared to a simple - multinomial logistic model, is clear when the delay between waves is not - identical for each individual. Also, if a individual missed an - intermediate interview, the information is lost, but taken into - account using an interpolation or extrapolation. - - hPijx is the probability to be observed in state i at age x+h - conditional to the observed state i at age x. The delay 'h' can be - split into an exact number (nh*stepm) of unobserved intermediate - states. This elementary transition (by month or quarter trimester, - semester or year) is model as a multinomial logistic. The hPx - matrix is simply the matrix product of nh*stepm elementary matrices - and the contribution of each individual to the likelihood is simply - hPijx. - - Also this programme outputs the covariance matrix of the parameters but also - of the life expectancies. It also computes the prevalence limits. - - Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). - Institut national d'études démographiques, Paris. - This software have been partly granted by Euro-REVES, a concerted action - from the European Union. - It is copyrighted identically to a GNU software product, ie programme and - software can be distributed freely for non commercial use. Latest version - can be accessed at http://euroreves.ined.fr/imach . - **********************************************************************/ - -#include -#include -#include -#include - -#define MAXLINE 256 -#define GNUPLOTPROGRAM "gnuplot" -/*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ -#define FILENAMELENGTH 80 -/*#define DEBUG*/ -#define windows -#define GLOCK_ERROR_NOPATH -1 /* empty path */ -#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ - -#define MAXPARM 30 /* Maximum number of parameters for the optimization */ -#define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */ - -#define NINTERVMAX 8 -#define NLSTATEMAX 8 /* Maximum number of live states (for func) */ -#define NDEATHMAX 8 /* Maximum number of dead states (for func) */ -#define NCOVMAX 8 /* Maximum number of covariates */ -#define MAXN 20000 -#define YEARM 12. /* Number of months per year */ -#define AGESUP 130 -#define AGEBASE 40 - - -int erreur; /* Error number */ -int nvar; -int cptcovn, cptcovage=0, cptcoveff=0,cptcov; -int npar=NPARMAX; -int nlstate=2; /* Number of live states */ -int ndeath=1; /* Number of dead states */ -int ncovmodel, ncovcol; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ -int popbased=0; - -int *wav; /* Number of waves for this individuual 0 is possible */ -int maxwav; /* Maxim number of waves */ -int jmin, jmax; /* min, max spacing between 2 waves */ -int mle, weightopt; -int **mw; /* mw[mi][i] is number of the mi wave for this individual */ -int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ -double jmean; /* Mean space between 2 waves */ -double **oldm, **newm, **savm; /* Working pointers to matrices */ -double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ -FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; -FILE *ficgp,*ficresprob,*ficpop; -FILE *ficreseij; - char filerese[FILENAMELENGTH]; - FILE *ficresvij; - char fileresv[FILENAMELENGTH]; - FILE *ficresvpl; - char fileresvpl[FILENAMELENGTH]; - -#define NR_END 1 -#define FREE_ARG char* -#define FTOL 1.0e-10 - -#define NRANSI -#define ITMAX 200 - -#define TOL 2.0e-4 - -#define CGOLD 0.3819660 -#define ZEPS 1.0e-10 -#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); - -#define GOLD 1.618034 -#define GLIMIT 100.0 -#define TINY 1.0e-20 - -static double maxarg1,maxarg2; -#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2)) -#define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2)) - -#define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a)) -#define rint(a) floor(a+0.5) - -static double sqrarg; -#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg) -#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} - -int imx; -int stepm; -/* Stepm, step in month: minimum step interpolation*/ - -int estepm; -/* Estepm, step in month to interpolate survival function in order to approximate Life Expectancy*/ - -int m,nb; -int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage; -double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; -double **pmmij, ***probs, ***mobaverage; -double dateintmean=0; - -double *weight; -int **s; /* Status */ -double *agedc, **covar, idx; -int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; - -double ftol=FTOL; /* Tolerance for computing Max Likelihood */ -double ftolhess; /* Tolerance for computing hessian */ - -/**************** split *************************/ -static int split( char *path, char *dirc, char *name, char *ext, char *finame ) -{ - char *s; /* pointer */ - int l1, l2; /* length counters */ - - l1 = strlen( path ); /* length of path */ - if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); -#ifdef windows - s = strrchr( path, '\\' ); /* find last / */ -#else - s = strrchr( path, '/' ); /* find last / */ -#endif - if ( s == NULL ) { /* no directory, so use current */ -#if defined(__bsd__) /* get current working directory */ - extern char *getwd( ); - - if ( getwd( dirc ) == NULL ) { -#else - extern char *getcwd( ); - - if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { -#endif - return( GLOCK_ERROR_GETCWD ); - } - strcpy( name, path ); /* we've got it */ - } else { /* strip direcotry from path */ - s++; /* after this, the filename */ - l2 = strlen( s ); /* length of filename */ - if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); - strcpy( name, s ); /* save file name */ - strncpy( dirc, path, l1 - l2 ); /* now the directory */ - dirc[l1-l2] = 0; /* add zero */ - } - l1 = strlen( dirc ); /* length of directory */ -#ifdef windows - if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; } -#else - if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; } -#endif - s = strrchr( name, '.' ); /* find last / */ - s++; - strcpy(ext,s); /* save extension */ - l1= strlen( name); - l2= strlen( s)+1; - strncpy( finame, name, l1-l2); - finame[l1-l2]= 0; - return( 0 ); /* we're done */ -} - - -/******************************************/ - -void replace(char *s, char*t) -{ - int i; - int lg=20; - i=0; - lg=strlen(t); - for(i=0; i<= lg; i++) { - (s[i] = t[i]); - if (t[i]== '\\') s[i]='/'; - } -} - -int nbocc(char *s, char occ) -{ - int i,j=0; - int lg=20; - i=0; - lg=strlen(s); - for(i=0; i<= lg; i++) { - if (s[i] == occ ) j++; - } - return j; -} - -void cutv(char *u,char *v, char*t, char occ) -{ - int i,lg,j,p=0; - i=0; - for(j=0; j<=strlen(t)-1; j++) { - if((t[j]!= occ) && (t[j+1]== occ)) p=j+1; - } - - lg=strlen(t); - for(j=0; j=(p+1))(v[j-p-1] = t[j]); - } -} - -/********************** nrerror ********************/ - -void nrerror(char error_text[]) -{ - fprintf(stderr,"ERREUR ...\n"); - fprintf(stderr,"%s\n",error_text); - exit(1); -} -/*********************** vector *******************/ -double *vector(int nl, int nh) -{ - double *v; - v=(double *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(double))); - if (!v) nrerror("allocation failure in vector"); - return v-nl+NR_END; -} - -/************************ free vector ******************/ -void free_vector(double*v, int nl, int nh) -{ - free((FREE_ARG)(v+nl-NR_END)); -} - -/************************ivector *******************************/ -int *ivector(long nl,long nh) -{ - int *v; - v=(int *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(int))); - if (!v) nrerror("allocation failure in ivector"); - return v-nl+NR_END; -} - -/******************free ivector **************************/ -void free_ivector(int *v, long nl, long nh) -{ - free((FREE_ARG)(v+nl-NR_END)); -} - -/******************* imatrix *******************************/ -int **imatrix(long nrl, long nrh, long ncl, long nch) - /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ -{ - long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; - int **m; - - /* allocate pointers to rows */ - m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*))); - if (!m) nrerror("allocation failure 1 in matrix()"); - m += NR_END; - m -= nrl; - - - /* allocate rows and set pointers to them */ - m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int))); - if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); - m[nrl] += NR_END; - m[nrl] -= ncl; - - for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; - - /* return pointer to array of pointers to rows */ - return m; -} - -/****************** free_imatrix *************************/ -void free_imatrix(m,nrl,nrh,ncl,nch) - int **m; - long nch,ncl,nrh,nrl; - /* free an int matrix allocated by imatrix() */ -{ - free((FREE_ARG) (m[nrl]+ncl-NR_END)); - free((FREE_ARG) (m+nrl-NR_END)); -} - -/******************* matrix *******************************/ -double **matrix(long nrl, long nrh, long ncl, long nch) -{ - long i, nrow=nrh-nrl+1, ncol=nch-ncl+1; - double **m; - - m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); - if (!m) nrerror("allocation failure 1 in matrix()"); - m += NR_END; - m -= nrl; - - m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); - if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); - m[nrl] += NR_END; - m[nrl] -= ncl; - - for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol; - return m; -} - -/*************************free matrix ************************/ -void free_matrix(double **m, long nrl, long nrh, long ncl, long nch) -{ - free((FREE_ARG)(m[nrl]+ncl-NR_END)); - free((FREE_ARG)(m+nrl-NR_END)); -} - -/******************* ma3x *******************************/ -double ***ma3x(long nrl, long nrh, long ncl, long nch, long nll, long nlh) -{ - long i, j, nrow=nrh-nrl+1, ncol=nch-ncl+1, nlay=nlh-nll+1; - double ***m; - - m=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double*))); - if (!m) nrerror("allocation failure 1 in matrix()"); - m += NR_END; - m -= nrl; - - m[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); - if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); - m[nrl] += NR_END; - m[nrl] -= ncl; - - for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol; - - m[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*nlay+NR_END)*sizeof(double))); - if (!m[nrl][ncl]) nrerror("allocation failure 3 in matrix()"); - m[nrl][ncl] += NR_END; - m[nrl][ncl] -= nll; - for (j=ncl+1; j<=nch; j++) - m[nrl][j]=m[nrl][j-1]+nlay; - - for (i=nrl+1; i<=nrh; i++) { - m[i][ncl]=m[i-1l][ncl]+ncol*nlay; - for (j=ncl+1; j<=nch; j++) - m[i][j]=m[i][j-1]+nlay; - } - return m; -} - -/*************************free ma3x ************************/ -void free_ma3x(double ***m, long nrl, long nrh, long ncl, long nch,long nll, long nlh) -{ - free((FREE_ARG)(m[nrl][ncl]+ nll-NR_END)); - free((FREE_ARG)(m[nrl]+ncl-NR_END)); - free((FREE_ARG)(m+nrl-NR_END)); -} - -/***************** f1dim *************************/ -extern int ncom; -extern double *pcom,*xicom; -extern double (*nrfunc)(double []); - -double f1dim(double x) -{ - int j; - double f; - double *xt; - - xt=vector(1,ncom); - for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j]; - f=(*nrfunc)(xt); - free_vector(xt,1,ncom); - return f; -} - -/*****************brent *************************/ -double brent(double ax, double bx, double cx, double (*f)(double), double tol, double *xmin) -{ - int iter; - double a,b,d,etemp; - double fu,fv,fw,fx; - double ftemp; - double p,q,r,tol1,tol2,u,v,w,x,xm; - double e=0.0; - - a=(ax < cx ? ax : cx); - b=(ax > cx ? ax : cx); - x=w=v=bx; - fw=fv=fx=(*f)(x); - for (iter=1;iter<=ITMAX;iter++) { - xm=0.5*(a+b); - tol2=2.0*(tol1=tol*fabs(x)+ZEPS); - /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/ - printf(".");fflush(stdout); -#ifdef DEBUG - printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol); - /* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */ -#endif - if (fabs(x-xm) <= (tol2-0.5*(b-a))){ - *xmin=x; - return fx; - } - ftemp=fu; - if (fabs(e) > tol1) { - r=(x-w)*(fx-fv); - q=(x-v)*(fx-fw); - p=(x-v)*q-(x-w)*r; - q=2.0*(q-r); - if (q > 0.0) p = -p; - q=fabs(q); - etemp=e; - e=d; - if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) - d=CGOLD*(e=(x >= xm ? a-x : b-x)); - else { - d=p/q; - u=x+d; - if (u-a < tol2 || b-u < tol2) - d=SIGN(tol1,xm-x); - } - } else { - d=CGOLD*(e=(x >= xm ? a-x : b-x)); - } - u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); - fu=(*f)(u); - if (fu <= fx) { - if (u >= x) a=x; else b=x; - SHFT(v,w,x,u) - SHFT(fv,fw,fx,fu) - } else { - if (u < x) a=u; else b=u; - if (fu <= fw || w == x) { - v=w; - w=u; - fv=fw; - fw=fu; - } else if (fu <= fv || v == x || v == w) { - v=u; - fv=fu; - } - } - } - nrerror("Too many iterations in brent"); - *xmin=x; - return fx; -} - -/****************** mnbrak ***********************/ - -void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc, - double (*func)(double)) -{ - double ulim,u,r,q, dum; - double fu; - - *fa=(*func)(*ax); - *fb=(*func)(*bx); - if (*fb > *fa) { - SHFT(dum,*ax,*bx,dum) - SHFT(dum,*fb,*fa,dum) - } - *cx=(*bx)+GOLD*(*bx-*ax); - *fc=(*func)(*cx); - while (*fb > *fc) { - r=(*bx-*ax)*(*fb-*fc); - q=(*bx-*cx)*(*fb-*fa); - u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ - (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); - ulim=(*bx)+GLIMIT*(*cx-*bx); - if ((*bx-u)*(u-*cx) > 0.0) { - fu=(*func)(u); - } else if ((*cx-u)*(u-ulim) > 0.0) { - fu=(*func)(u); - if (fu < *fc) { - SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) - SHFT(*fb,*fc,fu,(*func)(u)) - } - } else if ((u-ulim)*(ulim-*cx) >= 0.0) { - u=ulim; - fu=(*func)(u); - } else { - u=(*cx)+GOLD*(*cx-*bx); - fu=(*func)(u); - } - SHFT(*ax,*bx,*cx,u) - SHFT(*fa,*fb,*fc,fu) - } -} - -/*************** linmin ************************/ - -int ncom; -double *pcom,*xicom; -double (*nrfunc)(double []); - -void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [])) -{ - double brent(double ax, double bx, double cx, - double (*f)(double), double tol, double *xmin); - double f1dim(double x); - void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, - double *fc, double (*func)(double)); - int j; - double xx,xmin,bx,ax; - double fx,fb,fa; - - ncom=n; - pcom=vector(1,n); - xicom=vector(1,n); - nrfunc=func; - for (j=1;j<=n;j++) { - pcom[j]=p[j]; - xicom[j]=xi[j]; - } - ax=0.0; - xx=1.0; - mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); - *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); -#ifdef DEBUG - printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); -#endif - for (j=1;j<=n;j++) { - xi[j] *= xmin; - p[j] += xi[j]; - } - free_vector(xicom,1,n); - free_vector(pcom,1,n); -} - -/*************** powell ************************/ -void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, - double (*func)(double [])) -{ - void linmin(double p[], double xi[], int n, double *fret, - double (*func)(double [])); - int i,ibig,j; - double del,t,*pt,*ptt,*xit; - double fp,fptt; - double *xits; - pt=vector(1,n); - ptt=vector(1,n); - xit=vector(1,n); - xits=vector(1,n); - *fret=(*func)(p); - for (j=1;j<=n;j++) pt[j]=p[j]; - for (*iter=1;;++(*iter)) { - fp=(*fret); - ibig=0; - del=0.0; - printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret); - for (i=1;i<=n;i++) - printf(" %d %.12f",i, p[i]); - printf("\n"); - for (i=1;i<=n;i++) { - for (j=1;j<=n;j++) xit[j]=xi[j][i]; - fptt=(*fret); -#ifdef DEBUG - printf("fret=%lf \n",*fret); -#endif - printf("%d",i);fflush(stdout); - linmin(p,xit,n,fret,func); - if (fabs(fptt-(*fret)) > del) { - del=fabs(fptt-(*fret)); - ibig=i; - } -#ifdef DEBUG - printf("%d %.12e",i,(*fret)); - for (j=1;j<=n;j++) { - xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5); - printf(" x(%d)=%.12e",j,xit[j]); - } - for(j=1;j<=n;j++) - printf(" p=%.12e",p[j]); - printf("\n"); -#endif - } - if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { -#ifdef DEBUG - int k[2],l; - k[0]=1; - k[1]=-1; - printf("Max: %.12e",(*func)(p)); - for (j=1;j<=n;j++) - printf(" %.12e",p[j]); - printf("\n"); - for(l=0;l<=1;l++) { - for (j=1;j<=n;j++) { - ptt[j]=p[j]+(p[j]-pt[j])*k[l]; - printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]); - } - printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p))); - } -#endif - - - free_vector(xit,1,n); - free_vector(xits,1,n); - free_vector(ptt,1,n); - free_vector(pt,1,n); - return; - } - if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); - for (j=1;j<=n;j++) { - ptt[j]=2.0*p[j]-pt[j]; - xit[j]=p[j]-pt[j]; - pt[j]=p[j]; - } - fptt=(*func)(ptt); - if (fptt < fp) { - t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); - if (t < 0.0) { - linmin(p,xit,n,fret,func); - for (j=1;j<=n;j++) { - xi[j][ibig]=xi[j][n]; - xi[j][n]=xit[j]; - } -#ifdef DEBUG - printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); - for(j=1;j<=n;j++) - printf(" %.12e",xit[j]); - printf("\n"); -#endif - } - } - } -} - -/**** Prevalence limit ****************/ - -double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) -{ - /* Computes the prevalence limit in each live state at age x by left multiplying the unit - matrix by transitions matrix until convergence is reached */ - - int i, ii,j,k; - double min, max, maxmin, maxmax,sumnew=0.; - double **matprod2(); - double **out, cov[NCOVMAX], **pmij(); - double **newm; - double agefin, delaymax=50 ; /* Max number of years to converge */ - - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - } - - cov[1]=1.; - - /* Even if hstepm = 1, at least one multiplication by the unit matrix */ - for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){ - newm=savm; - /* Covariates have to be included here again */ - cov[2]=agefin; - - for (k=1; k<=cptcovn;k++) { - cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; - /* printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/ - } - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; - for (k=1; k<=cptcovprod;k++) - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; - - /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ - /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ - /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ - out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); - - savm=oldm; - oldm=newm; - maxmax=0.; - for(j=1;j<=nlstate;j++){ - min=1.; - max=0.; - for(i=1; i<=nlstate; i++) { - sumnew=0; - for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; - prlim[i][j]= newm[i][j]/(1-sumnew); - max=FMAX(max,prlim[i][j]); - min=FMIN(min,prlim[i][j]); - } - maxmin=max-min; - maxmax=FMAX(maxmax,maxmin); - } - if(maxmax < ftolpl){ - return prlim; - } - } -} - -/*************** transition probabilities ***************/ - -double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) -{ - double s1, s2; - /*double t34;*/ - int i,j,j1, nc, ii, jj; - - for(i=1; i<= nlstate; i++){ - for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/ - } - ps[i][j]=s2; - } - } - /*ps[3][2]=1;*/ - - for(i=1; i<= nlstate; i++){ - s1=0; - for(j=1; ji) { - printf(".%d%d",i,j);fflush(stdout); - hess[i][j]=hessij(p,delti,i,j); - hess[j][i]=hess[i][j]; - /*printf(" %lf ",hess[i][j]);*/ - } - } - } - printf("\n"); - - printf("\nInverting the hessian to get the covariance matrix. Wait...\n"); - - a=matrix(1,npar,1,npar); - y=matrix(1,npar,1,npar); - x=vector(1,npar); - indx=ivector(1,npar); - for (i=1;i<=npar;i++) - for (j=1;j<=npar;j++) a[i][j]=hess[i][j]; - ludcmp(a,npar,indx,&pd); - - for (j=1;j<=npar;j++) { - for (i=1;i<=npar;i++) x[i]=0; - x[j]=1; - lubksb(a,npar,indx,x); - for (i=1;i<=npar;i++){ - matcov[i][j]=x[i]; - } - } - - printf("\n#Hessian matrix#\n"); - for (i=1;i<=npar;i++) { - for (j=1;j<=npar;j++) { - printf("%.3e ",hess[i][j]); - } - printf("\n"); - } - - /* Recompute Inverse */ - for (i=1;i<=npar;i++) - for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; - ludcmp(a,npar,indx,&pd); - - /* printf("\n#Hessian matrix recomputed#\n"); - - for (j=1;j<=npar;j++) { - for (i=1;i<=npar;i++) x[i]=0; - x[j]=1; - lubksb(a,npar,indx,x); - for (i=1;i<=npar;i++){ - y[i][j]=x[i]; - printf("%.3e ",y[i][j]); - } - printf("\n"); - } - */ - - free_matrix(a,1,npar,1,npar); - free_matrix(y,1,npar,1,npar); - free_vector(x,1,npar); - free_ivector(indx,1,npar); - free_matrix(hess,1,npar,1,npar); - - -} - -/*************** hessian matrix ****************/ -double hessii( double x[], double delta, int theta, double delti[]) -{ - int i; - int l=1, lmax=20; - double k1,k2; - double p2[NPARMAX+1]; - double res; - double delt, delts, nkhi=10.,nkhif=1., khi=1.e-4; - double fx; - int k=0,kmax=10; - double l1; - - fx=func(x); - for (i=1;i<=npar;i++) p2[i]=x[i]; - for(l=0 ; l <=lmax; l++){ - l1=pow(10,l); - delts=delt; - for(k=1 ; k khi/nkhif) || (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */ - k=kmax; l=lmax*10.; - } - else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ - delts=delt; - } - } - } - delti[theta]=delts; - return res; - -} - -double hessij( double x[], double delti[], int thetai,int thetaj) -{ - int i; - int l=1, l1, lmax=20; - double k1,k2,k3,k4,res,fx; - double p2[NPARMAX+1]; - int k; - - fx=func(x); - for (k=1; k<=2; k++) { - for (i=1;i<=npar;i++) p2[i]=x[i]; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; - k1=func(p2)-fx; - - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; - k2=func(p2)-fx; - - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; - k3=func(p2)-fx; - - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; - k4=func(p2)-fx; - res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */ -#ifdef DEBUG - 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); -#endif - } - return res; -} - -/************** Inverse of matrix **************/ -void ludcmp(double **a, int n, int *indx, double *d) -{ - int i,imax,j,k; - double big,dum,sum,temp; - double *vv; - - vv=vector(1,n); - *d=1.0; - for (i=1;i<=n;i++) { - 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"); - vv[i]=1.0/big; - } - for (j=1;j<=n;j++) { - for (i=1;i= big) { - big=dum; - imax=i; - } - } - if (j != imax) { - for (k=1;k<=n;k++) { - dum=a[imax][k]; - a[imax][k]=a[j][k]; - a[j][k]=dum; - } - *d = -(*d); - vv[imax]=vv[j]; - } - indx[j]=imax; - if (a[j][j] == 0.0) a[j][j]=TINY; - if (j != n) { - dum=1.0/(a[j][j]); - for (i=j+1;i<=n;i++) a[i][j] *= dum; - } - } - free_vector(vv,1,n); /* Doesn't work */ -; -} - -void lubksb(double **a, int n, int *indx, double b[]) -{ - int i,ii=0,ip,j; - double sum; - - for (i=1;i<=n;i++) { - ip=indx[i]; - sum=b[ip]; - b[ip]=b[i]; - if (ii) - for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; - else if (sum) ii=i; - b[i]=sum; - } - for (i=n;i>=1;i--) { - sum=b[i]; - for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; - b[i]=sum/a[i][i]; - } -} - -/************ Frequencies ********************/ -void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2) -{ /* Some frequencies */ - - int i, m, jk, k1,i1, j1, bool, z1,z2,j; - double ***freq; /* Frequencies */ - double *pp; - double pos, k2, dateintsum=0,k2cpt=0; - FILE *ficresp; - char fileresp[FILENAMELENGTH]; - - pp=vector(1,nlstate); - probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); - strcpy(fileresp,"p"); - strcat(fileresp,fileres); - if((ficresp=fopen(fileresp,"w"))==NULL) { - printf("Problem with prevalence resultfile: %s\n", fileresp); - exit(0); - } - freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); - j1=0; - - j=cptcoveff; - if (cptcovn<1) {j=1;ncodemax[1]=1;} - - for(k1=1; k1<=j;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ - j1++; - /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); - scanf("%d", i);*/ - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) - for(m=agemin; m <= agemax+3; m++) - freq[i][jk][m]=0; - - dateintsum=0; - k2cpt=0; - for (i=1; i<=imx; i++) { - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) - bool=0; - } - if (bool==1) { - for(m=firstpass; m<=lastpass; m++){ - k2=anint[m][i]+(mint[m][i]/12.); - if ((k2>=dateprev1) && (k2<=dateprev2)) { - if(agev[m][i]==0) agev[m][i]=agemax+1; - if(agev[m][i]==1) agev[m][i]=agemax+2; - if (m1) && (agev[m][i]< (agemax+3))) { - dateintsum=dateintsum+k2; - k2cpt++; - } - } - } - } - } - - fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); - - if (cptcovn>0) { - fprintf(ficresp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficresp, "**********\n#"); - } - for(i=1; i<=nlstate;i++) - fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); - fprintf(ficresp, "\n"); - - for(i=(int)agemin; i <= (int)agemax+3; i++){ - if(i==(int)agemax+3) - printf("Total"); - else - printf("Age %d", i); - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][i]; - if(pp[jk]>=1.e-10) - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); - else - printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - } - - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - pp[jk] += freq[jk][m][i]; - } - - for(jk=1,pos=0; jk <=nlstate ; jk++) - pos += pp[jk]; - for(jk=1; jk <=nlstate ; jk++){ - if(pos>=1.e-5) - printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - else - printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - if( i <= (int) agemax){ - if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, 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]);*/ - } - else - fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos); - } - } - - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) - if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); - if(i <= (int) agemax) - fprintf(ficresp,"\n"); - printf("\n"); - } - } - } - dateintmean=dateintsum/k2cpt; - - fclose(ficresp); - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); - free_vector(pp,1,nlstate); - - /* End of Freq */ -} - -/************ Prevalence ********************/ -void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate) -{ /* Some frequencies */ - - int i, m, jk, k1, i1, j1, bool, z1,z2,j; - double ***freq; /* Frequencies */ - double *pp; - double pos, k2; - - pp=vector(1,nlstate); - probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); - - freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); - j1=0; - - j=cptcoveff; - if (cptcovn<1) {j=1;ncodemax[1]=1;} - - for(k1=1; k1<=j;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ - j1++; - - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) - for(m=agemin; m <= agemax+3; m++) - freq[i][jk][m]=0; - - for (i=1; i<=imx; i++) { - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) - bool=0; - } - if (bool==1) { - for(m=firstpass; m<=lastpass; m++){ - k2=anint[m][i]+(mint[m][i]/12.); - if ((k2>=dateprev1) && (k2<=dateprev2)) { - if(agev[m][i]==0) agev[m][i]=agemax+1; - if(agev[m][i]==1) agev[m][i]=agemax+2; - if (m0) - freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i]; - else - freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; - freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; - } - } - } - } - } - for(i=(int)agemin; i <= (int)agemax+3; i++){ - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][i]; - } - - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - pp[jk] += freq[jk][m][i]; - } - - for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk]; - - for(jk=1; jk <=nlstate ; jk++){ - if( i <= (int) agemax){ - if(pos>=1.e-5){ - probs[i][jk][j1]= pp[jk]/pos; - } - } - } - - } - } - } - - - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); - free_vector(pp,1,nlstate); - -} /* End of Freq */ - -/************* Waves Concatenation ***************/ - -void concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm) -{ - /* 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 mi (mi=1 to wav[i]) effective wave of individual i - dh[m][i] of 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. - */ - - int i, mi, m; - /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; - double sum=0., jmean=0.;*/ - - int j, k=0,jk, ju, jl; - double sum=0.; - jmin=1e+5; - jmax=-1; - jmean=0.; - for(i=1; i<=imx; i++){ - mi=0; - m=firstpass; - while(s[m][i] <= nlstate){ - if(s[m][i]>=1) - mw[++mi][i]=m; - if(m >=lastpass) - break; - else - m++; - }/* end while */ - if (s[m][i] > nlstate){ - mi++; /* Death is another wave */ - /* if(mi==0) never been interviewed correctly before death */ - /* Only death is a correct wave */ - mw[mi][i]=m; - } - - wav[i]=mi; - if(mi==0) - printf("Warning, no any valid information for:%d line=%d\n",num[i],i); - } - - for(i=1; i<=imx; i++){ - for(mi=1; mi nlstate) { - 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 */ - k=k+1; - if (j >= jmax) jmax=j; - if (j <= jmin) jmin=j; - sum=sum+j; - /*if (j<0) printf("j=%d num=%d \n",j,i); */ - } - } - else{ - j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); - k=k+1; - if (j >= jmax) jmax=j; - else if (j <= jmin)jmin=j; - /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ - sum=sum+j; - } - jk= j/stepm; - jl= j -jk*stepm; - ju= j -(jk+1)*stepm; - if(jl <= -ju) - dh[mi][i]=jk; - else - dh[mi][i]=jk+1; - if(dh[mi][i]==0) - dh[mi][i]=1; /* At least one step */ - } - } - } - jmean=sum/k; - printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); - } -/*********** Tricode ****************************/ -void tricode(int *Tvar, int **nbcode, int imx) -{ - int Ndum[20],ij=1, k, j, i; - int cptcode=0; - cptcoveff=0; - - for (k=0; k<19; k++) Ndum[k]=0; - for (k=1; k<=7; k++) ncodemax[k]=0; - - for (j=1; j<=(cptcovn+2*cptcovprod); j++) { - for (i=1; i<=imx; i++) { - ij=(int)(covar[Tvar[j]][i]); - Ndum[ij]++; - /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ - if (ij > cptcode) cptcode=ij; - } - - for (i=0; i<=cptcode; i++) { - if(Ndum[i]!=0) ncodemax[j]++; - } - ij=1; - - - for (i=1; i<=ncodemax[j]; i++) { - for (k=0; k<=19; k++) { - if (Ndum[k] != 0) { - nbcode[Tvar[j]][ij]=k; - - ij++; - } - if (ij > ncodemax[j]) break; - } - } - } - - for (k=0; k<19; k++) Ndum[k]=0; - - for (i=1; i<=ncovmodel-2; i++) { - ij=Tvar[i]; - Ndum[ij]++; - } - - ij=1; - for (i=1; i<=10; i++) { - if((Ndum[i]!=0) && (i<=ncovcol)){ - Tvaraff[ij]=i; - ij++; - } - } - - cptcoveff=ij-1; -} - -/*********** Health Expectancies ****************/ - -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov ) - -{ - /* Health expectancies */ - int i, j, nhstepm, hstepm, h, nstepm, k, cptj; - double age, agelim, hf; - double ***p3mat,***varhe; - double **dnewm,**doldm; - double *xp; - double **gp, **gm; - double ***gradg, ***trgradg; - int theta; - - varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage); - xp=vector(1,npar); - dnewm=matrix(1,nlstate*2,1,npar); - doldm=matrix(1,nlstate*2,1,nlstate*2); - - fprintf(ficreseij,"# Health expectancies\n"); - fprintf(ficreseij,"# Age"); - for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate;j++) - fprintf(ficreseij," %1d-%1d (SE)",i,j); - fprintf(ficreseij,"\n"); - - if(estepm < stepm){ - printf ("Problem %d lower than %d\n",estepm, stepm); - } - else hstepm=estepm; - /* We compute the life expectancy from trapezoids spaced every estepm months - * This is mainly to measure the difference between two models: for example - * if stepm=24 months pijx are given only every 2 years and by summing them - * we are calculating an estimate of the Life Expectancy assuming a linear - * progression inbetween and thus overestimating or underestimating according - * to the curvature of the survival function. If, for the same date, we - * estimate the model with stepm=1 month, we can keep estepm to 24 months - * to compare the new estimate of Life expectancy with the same linear - * hypothesis. A more precise result, taking into account a more precise - * curvature will be obtained if estepm is as small as stepm. */ - - /* 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. - nhstepm is the number of hstepm from age to agelim - nstepm is the number of stepm from age to agelin. - Look at hpijx to understand the reason of that which relies in memory size - and note for a fixed period like estepm months */ - /* 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 - means that if the survival funtion is printed only each two years of age and if - you sum them up and add 1 year (area under the trapezoids) you won't get the same - results. So we changed our mind and took the option of the best precision. - */ - hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ - - agelim=AGESUP; - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - /* nhstepm age range expressed in number of stepm */ - nstepm=(int) rint((agelim-age)*YEARM/stepm); - /* Typically if 20 years nstepm = 20*12/6=40 stepm */ - /* if (stepm >= YEARM) hstepm=1;*/ - nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2); - gp=matrix(0,nhstepm,1,nlstate*2); - gm=matrix(0,nhstepm,1,nlstate*2); - - /* Computed by stepm unit matrices, product of hstepm matrices, stored - in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij); - - - hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ - - /* Computing Variances of health expectancies */ - - for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++){ - xp[i] = x[i] + (i==theta ?delti[theta]:0); - } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - - cptj=0; - for(j=1; j<= nlstate; j++){ - for(i=1; i<=nlstate; i++){ - cptj=cptj+1; - for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){ - gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; - } - } - } - - - for(i=1; i<=npar; i++) - xp[i] = x[i] - (i==theta ?delti[theta]:0); - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - - cptj=0; - for(j=1; j<= nlstate; j++){ - for(i=1;i<=nlstate;i++){ - cptj=cptj+1; - for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){ - gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; - } - } - } - - - - for(j=1; j<= nlstate*2; j++) - for(h=0; h<=nhstepm-1; h++){ - gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; - } - - } - -/* End theta */ - - trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar); - - for(h=0; h<=nhstepm-1; h++) - for(j=1; j<=nlstate*2;j++) - for(theta=1; theta <=npar; theta++) - trgradg[h][j][theta]=gradg[h][theta][j]; - - - for(i=1;i<=nlstate*2;i++) - for(j=1;j<=nlstate*2;j++) - varhe[i][j][(int)age] =0.; - - printf("%d||",(int)age);fflush(stdout); - for(h=0;h<=nhstepm-1;h++){ - for(k=0;k<=nhstepm-1;k++){ - matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]); - for(i=1;i<=nlstate*2;i++) - for(j=1;j<=nlstate*2;j++) - varhe[i][j][(int)age] += doldm[i][j]*hf*hf; - } - } - - - /* Computing expectancies */ - for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate;j++) - for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ - eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; - -/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ - - } - - fprintf(ficreseij,"%3.0f",age ); - cptj=0; - for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate;j++){ - cptj++; - fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) ); - } - fprintf(ficreseij,"\n"); - - free_matrix(gm,0,nhstepm,1,nlstate*2); - free_matrix(gp,0,nhstepm,1,nlstate*2); - free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2); - free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar); - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - } - free_vector(xp,1,npar); - free_matrix(dnewm,1,nlstate*2,1,npar); - free_matrix(doldm,1,nlstate*2,1,nlstate*2); - free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage); -} - -/************ Variance ******************/ -void varevsij(char fileres[], 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 ij, int estepm) -{ - /* Variance of health expectancies */ - /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ - double **newm; - double **dnewm,**doldm; - int i, j, nhstepm, hstepm, h, nstepm ; - int k, cptcode; - double *xp; - double **gp, **gm; - double ***gradg, ***trgradg; - double ***p3mat; - double age,agelim, hf; - int theta; - - fprintf(ficresvij,"# Covariances of life expectancies\n"); - fprintf(ficresvij,"# Age"); - for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate;j++) - fprintf(ficresvij," Cov(e%1d, e%1d)",i,j); - fprintf(ficresvij,"\n"); - - xp=vector(1,npar); - dnewm=matrix(1,nlstate,1,npar); - doldm=matrix(1,nlstate,1,nlstate); - - if(estepm < stepm){ - printf ("Problem %d lower than %d\n",estepm, stepm); - } - else hstepm=estepm; - /* 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. - nhstepm is the number of hstepm from age to agelim - nstepm is the number of stepm from age to agelin. - Look at hpijx to understand the reason of that which relies in memory size - and note for a fixed period like k years */ - /* 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 - means that if the survival funtion is printed only each two years of age and if - you sum them up and add 1 year (area under the trapezoids) you won't get the same - results. So we changed our mind and took the option of the best precision. - */ - hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ - agelim = AGESUP; - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ - nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - gradg=ma3x(0,nhstepm,1,npar,1,nlstate); - gp=matrix(0,nhstepm,1,nlstate); - gm=matrix(0,nhstepm,1,nlstate); - - for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++){ /* Computes gradient */ - 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,ij); - - if (popbased==1) { - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][ij]; - } - - for(j=1; j<= nlstate; j++){ - for(h=0; h<=nhstepm; h++){ - for(i=1, gp[h][j]=0.;i<=nlstate;i++) - gp[h][j] += prlim[i][i]*p3mat[i][j][h]; - } - } - - for(i=1; i<=npar; i++) /* Computes gradient */ - 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,ij); - - if (popbased==1) { - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][ij]; - } - - for(j=1; j<= nlstate; j++){ - for(h=0; h<=nhstepm; h++){ - for(i=1, gm[h][j]=0.;i<=nlstate;i++) - gm[h][j] += prlim[i][i]*p3mat[i][j][h]; - } - } - - for(j=1; j<= nlstate; j++) - for(h=0; h<=nhstepm; h++){ - gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; - } - } /* End theta */ - - trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); - - for(h=0; h<=nhstepm; h++) - for(j=1; j<=nlstate;j++) - for(theta=1; theta <=npar; theta++) - trgradg[h][j][theta]=gradg[h][theta][j]; - - hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ - for(i=1;i<=nlstate;i++) - for(j=1;j<=nlstate;j++) - vareij[i][j][(int)age] =0.; - - for(h=0;h<=nhstepm;h++){ - for(k=0;k<=nhstepm;k++){ - matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]); - for(i=1;i<=nlstate;i++) - for(j=1;j<=nlstate;j++) - vareij[i][j][(int)age] += doldm[i][j]*hf*hf; - } - } - - fprintf(ficresvij,"%.0f ",age ); - for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate;j++){ - fprintf(ficresvij," %.4f", vareij[i][j][(int)age]); - } - fprintf(ficresvij,"\n"); - free_matrix(gp,0,nhstepm,1,nlstate); - free_matrix(gm,0,nhstepm,1,nlstate); - free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate); - free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar); - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - } /* End age */ - - free_vector(xp,1,npar); - free_matrix(doldm,1,nlstate,1,npar); - free_matrix(dnewm,1,nlstate,1,nlstate); - -} - -/************ 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 ij) -{ - /* Variance of prevalence limit */ - /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ - double **newm; - double **dnewm,**doldm; - int i, j, nhstepm, hstepm; - int k, cptcode; - double *xp; - double *gp, *gm; - double **gradg, **trgradg; - double age,agelim; - int theta; - - fprintf(ficresvpl,"# Standard deviation of prevalences limit\n"); - fprintf(ficresvpl,"# Age"); - for(i=1; i<=nlstate;i++) - fprintf(ficresvpl," %1d-%1d",i,i); - fprintf(ficresvpl,"\n"); - - xp=vector(1,npar); - dnewm=matrix(1,nlstate,1,npar); - doldm=matrix(1,nlstate,1,nlstate); - - hstepm=1*YEARM; /* Every year of age */ - hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */ - agelim = AGESUP; - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ - if (stepm >= YEARM) hstepm=1; - nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ - gradg=matrix(1,npar,1,nlstate); - gp=vector(1,nlstate); - gm=vector(1,nlstate); - - for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++){ /* Computes gradient */ - xp[i] = x[i] + (i==theta ?delti[theta]:0); - } - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - for(i=1;i<=nlstate;i++) - gp[i] = prlim[i][i]; - - for(i=1; i<=npar; i++) /* Computes gradient */ - xp[i] = x[i] - (i==theta ?delti[theta]:0); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - for(i=1;i<=nlstate;i++) - gm[i] = prlim[i][i]; - - for(i=1;i<=nlstate;i++) - gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; - } /* End theta */ - - trgradg =matrix(1,nlstate,1,npar); - - for(j=1; j<=nlstate;j++) - for(theta=1; theta <=npar; theta++) - trgradg[j][theta]=gradg[theta][j]; - - for(i=1;i<=nlstate;i++) - varpl[i][(int)age] =0.; - matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); - for(i=1;i<=nlstate;i++) - varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ - - fprintf(ficresvpl,"%.0f ",age ); - for(i=1; i<=nlstate;i++) - fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); - fprintf(ficresvpl,"\n"); - free_vector(gp,1,nlstate); - free_vector(gm,1,nlstate); - free_matrix(gradg,1,npar,1,nlstate); - free_matrix(trgradg,1,nlstate,1,npar); - } /* End age */ - - free_vector(xp,1,npar); - free_matrix(doldm,1,nlstate,1,npar); - free_matrix(dnewm,1,nlstate,1,nlstate); - -} - -/************ Variance of one-step probabilities ******************/ -void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax) -{ - int i, j, i1, k1, j1, z1; - int k=0, cptcode; - double **dnewm,**doldm; - double *xp; - double *gp, *gm; - double **gradg, **trgradg; - double age,agelim, cov[NCOVMAX]; - int theta; - char fileresprob[FILENAMELENGTH]; - - strcpy(fileresprob,"prob"); - strcat(fileresprob,fileres); - if((ficresprob=fopen(fileresprob,"w"))==NULL) { - printf("Problem with resultfile: %s\n", fileresprob); - } - printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); - -fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n"); - fprintf(ficresprob,"# Age"); - for(i=1; i<=nlstate;i++) - for(j=1; j<=(nlstate+ndeath);j++) - fprintf(ficresprob," p%1d-%1d (SE)",i,j); - - - fprintf(ficresprob,"\n"); - - - xp=vector(1,npar); - dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); - doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); - - cov[1]=1; - j=cptcoveff; - if (cptcovn<1) {j=1;ncodemax[1]=1;} - j1=0; - for(k1=1; k1<=1;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ - j1++; - - if (cptcovn>0) { - fprintf(ficresprob, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficresprob, "**********\n#"); - } - - for (age=bage; age<=fage; age ++){ - cov[2]=age; - for (k=1; k<=cptcovn;k++) { - cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]]; - - } - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; - for (k=1; k<=cptcovprod;k++) - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; - - gradg=matrix(1,npar,1,9); - trgradg=matrix(1,9,1,npar); - gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); - gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); - - for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++) - xp[i] = x[i] + (i==theta ?delti[theta]:0); - - pmij(pmmij,cov,ncovmodel,xp,nlstate); - - k=0; - for(i=1; i<= (nlstate+ndeath); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gp[k]=pmmij[i][j]; - } - } - - for(i=1; i<=npar; i++) - xp[i] = x[i] - (i==theta ?delti[theta]:0); - - pmij(pmmij,cov,ncovmodel,xp,nlstate); - k=0; - for(i=1; i<=(nlstate+ndeath); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gm[k]=pmmij[i][j]; - } - } - - for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) - gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; - } - - for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++) - for(theta=1; theta <=npar; theta++) - trgradg[j][theta]=gradg[theta][j]; - - matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg); - - pmij(pmmij,cov,ncovmodel,x,nlstate); - - k=0; - for(i=1; i<=(nlstate+ndeath); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gm[k]=pmmij[i][j]; - } - } - - /*printf("\n%d ",(int)age); - for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){ - printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); - }*/ - - fprintf(ficresprob,"\n%d ",(int)age); - - for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++) - fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i])); - - } - } - free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); - free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); - free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); - free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); - } - free_vector(xp,1,npar); - fclose(ficresprob); - -} - -/******************* Printing html file ***********/ -void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \ - int lastpass, int stepm, int weightopt, char model[],\ - int imx,int jmin, int jmax, double jmeanint,char optionfile[], \ - char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\ - char version[], int popforecast, int estepm ){ - int jj1, k1, i1, cpt; - FILE *fichtm; - /*char optionfilehtm[FILENAMELENGTH];*/ - - strcpy(optionfilehtm,optionfile); - strcat(optionfilehtm,".htm"); - if((fichtm=fopen(optionfilehtm,"w"))==NULL) { - printf("Problem with %s \n",optionfilehtm), exit(0); - } - - fprintf(fichtm," %s
\n -Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n -\n -Total number of observations=%d
\n -Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n -
-