model. More health states you consider, more time is necessary to reach the\r
Maximum Likelihood of the parameters involved in the model. The\r
simplest model is the multinomial logistic model where pij is the\r
- probabibility to be observed in state j at the second wave\r
+ probability to be observed in state j at the second wave\r
conditional to be observed in state i at the first wave. Therefore\r
the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where\r
'age' is age and 'sex' is a covariate. If you want to have a more\r
for (k=0; k<=19; k++) {\r
if (Ndum[k] != 0) {\r
nbcode[Tvar[j]][ij]=k; \r
- /* printf("nbcodeaaaaaaaaaaa=%d Tvar[j]=%d ij=%d j=%d",nbcode[Tvar[j]][ij],Tvar[j],ij,j);*/\r
+ \r
ij++;\r
}\r
if (ij > ncodemax[j]) break; \r
/* Computed by stepm unit matrices, product of hstepm matrices, stored\r
in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */\r
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij); \r
+ \r
+ /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/\r
+\r
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */\r
for(i=1; i<=nlstate;i++)\r
for(j=1; j<=nlstate;j++)\r
}\r
\r
/************ Variance of one-step probabilities ******************/\r
-void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)\r
+void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)\r
{\r
- int i, j;\r
+ int i, j, i1, k1, j1, z1;\r
int k=0, cptcode;\r
double **dnewm,**doldm;\r
double *xp;\r
doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));\r
\r
cov[1]=1;\r
- for (age=bage; age<=fage; age ++){ \r
- cov[2]=age;\r
- gradg=matrix(1,npar,1,9);\r
- trgradg=matrix(1,9,1,npar);\r
- gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));\r
- gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));\r
+ j=cptcoveff;\r
+ if (cptcovn<1) {j=1;ncodemax[1]=1;}\r
+ j1=0;\r
+ for(k1=1; k1<=1;k1++){\r
+ for(i1=1; i1<=ncodemax[k1];i1++){ \r
+ j1++;\r
+\r
+ if (cptcovn>0) {\r
+ fprintf(ficresprob, "\n#********** Variable "); \r
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);\r
+ fprintf(ficresprob, "**********\n#");\r
+ }\r
\r
- for(theta=1; theta <=npar; theta++){\r
- for(i=1; i<=npar; i++)\r
- xp[i] = x[i] + (i==theta ?delti[theta]:0);\r
- \r
- pmij(pmmij,cov,ncovmodel,xp,nlstate);\r
- \r
- k=0;\r
- for(i=1; i<= (nlstate+ndeath); i++){\r
- for(j=1; j<=(nlstate+ndeath);j++){\r
- k=k+1;\r
- gp[k]=pmmij[i][j];\r
+ for (age=bage; age<=fage; age ++){ \r
+ cov[2]=age;\r
+ for (k=1; k<=cptcovn;k++) {\r
+ cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];\r
+ \r
}\r
- }\r
-\r
- for(i=1; i<=npar; i++)\r
- xp[i] = x[i] - (i==theta ?delti[theta]:0);\r
+ for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];\r
+ for (k=1; k<=cptcovprod;k++)\r
+ cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];\r
+ \r
+ gradg=matrix(1,npar,1,9);\r
+ trgradg=matrix(1,9,1,npar);\r
+ gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));\r
+ gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));\r
\r
-\r
- pmij(pmmij,cov,ncovmodel,xp,nlstate);\r
- k=0;\r
- for(i=1; i<=(nlstate+ndeath); i++){\r
- for(j=1; j<=(nlstate+ndeath);j++){\r
- k=k+1;\r
- gm[k]=pmmij[i][j];\r
- }\r
- }\r
+ for(theta=1; theta <=npar; theta++){\r
+ for(i=1; i<=npar; i++)\r
+ xp[i] = x[i] + (i==theta ?delti[theta]:0);\r
+ \r
+ pmij(pmmij,cov,ncovmodel,xp,nlstate);\r
+ \r
+ k=0;\r
+ for(i=1; i<= (nlstate+ndeath); i++){\r
+ for(j=1; j<=(nlstate+ndeath);j++){\r
+ k=k+1;\r
+ gp[k]=pmmij[i][j];\r
+ }\r
+ }\r
+ \r
+ for(i=1; i<=npar; i++)\r
+ xp[i] = x[i] - (i==theta ?delti[theta]:0);\r
+ \r
+ pmij(pmmij,cov,ncovmodel,xp,nlstate);\r
+ k=0;\r
+ for(i=1; i<=(nlstate+ndeath); i++){\r
+ for(j=1; j<=(nlstate+ndeath);j++){\r
+ k=k+1;\r
+ gm[k]=pmmij[i][j];\r
+ }\r
+ }\r
\r
- for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) \r
- gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; \r
- }\r
-\r
- for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)\r
- for(theta=1; theta <=npar; theta++)\r
- trgradg[j][theta]=gradg[theta][j];\r
- \r
- matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);\r
- matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);\r
-\r
- pmij(pmmij,cov,ncovmodel,x,nlstate);\r
+ for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) \r
+ gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; \r
+ }\r
\r
- k=0;\r
- for(i=1; i<=(nlstate+ndeath); i++){\r
- for(j=1; j<=(nlstate+ndeath);j++){\r
- k=k+1;\r
- gm[k]=pmmij[i][j];\r
+ for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)\r
+ for(theta=1; theta <=npar; theta++)\r
+ trgradg[j][theta]=gradg[theta][j];\r
+ \r
+ matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);\r
+ matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);\r
+ \r
+ pmij(pmmij,cov,ncovmodel,x,nlstate);\r
+ \r
+ k=0;\r
+ for(i=1; i<=(nlstate+ndeath); i++){\r
+ for(j=1; j<=(nlstate+ndeath);j++){\r
+ k=k+1;\r
+ gm[k]=pmmij[i][j];\r
+ }\r
}\r
- }\r
\r
/*printf("\n%d ",(int)age);\r
for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){\r
- \r
-\r
printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));\r
}*/\r
\r
- fprintf(ficresprob,"\n%d ",(int)age);\r
-\r
- for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){\r
- if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);\r
-if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);\r
- }\r
+ fprintf(ficresprob,"\n%d ",(int)age);\r
\r
+ for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)\r
+ fprintf(ficresprob,"%.3e (%.3e) ",gm[i],doldm[i][i]);\r
+ \r
+ }\r
+ }\r
free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));\r
free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));\r
free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);\r
free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);\r
-}\r
- free_vector(xp,1,npar);\r
-fclose(ficresprob);\r
-\r
+ }\r
+ free_vector(xp,1,npar);\r
+ fclose(ficresprob);\r
+ \r
}\r
\r
/******************* Printing html file ***********/\r
if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;\r
if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;\r
}*/\r
- \r
- /* for (i=1; i<=imx; i++){\r
+ /* for (i=1; i<=imx; i++){\r
if (s[4][i]==9) s[4][i]=-1; \r
- 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]));}\r
- */\r
+ 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]));}*/\r
+ \r
\r
/* Calculation of the number of parameter from char model*/\r
Tvar=ivector(1,15); \r
}\r
}\r
\r
- /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/\r
+ varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);\r
\r
fclose(ficrespij);\r
\r