From 1b1e97b93f674d872be751cb438712150d5ddcf4 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Agn=C3=A8s=20Li=C3=A8vre?= Date: Fri, 5 Apr 2002 15:45:00 +0000 Subject: [PATCH] *** empty log message *** --- src/imach.c | 158 +++++++++++++++++++++++++++++----------------------- 1 file changed, 89 insertions(+), 69 deletions(-) diff --git a/src/imach.c b/src/imach.c index 28c1ec7..5ba0474 100644 --- a/src/imach.c +++ b/src/imach.c @@ -14,7 +14,7 @@ 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 - probabibility to be observed in state j at the second wave + 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 @@ -1503,7 +1503,7 @@ void tricode(int *Tvar, int **nbcode, int imx) for (k=0; k<=19; k++) { if (Ndum[k] != 0) { nbcode[Tvar[j]][ij]=k; - /* printf("nbcodeaaaaaaaaaaa=%d Tvar[j]=%d ij=%d j=%d",nbcode[Tvar[j]][ij],Tvar[j],ij,j);*/ + ij++; } if (ij > ncodemax[j]) break; @@ -1585,6 +1585,9 @@ void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, in /* 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); + + /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/ + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) @@ -1822,9 +1825,9 @@ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], dou } /************ Variance of one-step probabilities ******************/ -void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij) +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; + int i, j, i1, k1, j1, z1; int k=0, cptcode; double **dnewm,**doldm; double *xp; @@ -1847,83 +1850,101 @@ void varprob(char fileres[], double **matcov, double x[], double delti[], int nl doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); cov[1]=1; - for (age=bage; age<=fage; age ++){ - cov[2]=age; - 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)); + 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(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 (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(i=1; i<=npar; i++) - xp[i] = x[i] - (i==theta ?delti[theta]:0); + 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)); - - 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(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); + for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) + gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; + } - 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(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++){ - if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]); -if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],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],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); - + } + free_vector(xp,1,npar); + fclose(ficresprob); + } /******************* Printing html file ***********/ @@ -2751,11 +2772,10 @@ while((c=getc(ficpar))=='#' && c!= EOF){ if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3; if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3; }*/ - - /* for (i=1; i<=imx; i++){ + /* for (i=1; i<=imx; i++){ if (s[4][i]==9) s[4][i]=-1; - printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));} - */ + printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/ + /* Calculation of the number of parameter from char model*/ Tvar=ivector(1,15); @@ -3220,7 +3240,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } } - /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/ + varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax); fclose(ficrespij); -- 2.43.0