/* $Id$
$State$
$Log$
+ Revision 1.318 2022/05/24 08:10:59 brouard
+ * imach.c (Module): Some attempts to find a bug of wrong estimates
+ of confidencce intervals with product in the equation modelC
+
Revision 1.317 2022/05/15 15:06:23 brouard
* imach.c (Module): Some minor improvements
#define POWELLNOF3INFF1TEST /* Skip test */
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
+/* #define FLATSUP *//* Suppresses directions where likelihood is flat */
#include <math.h>
#include <stdio.h>
double ***cotqvar; /* Time varying quantitative covariate itqv */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+/* Some documentation */
+ /* Design original data
+ * V1 V2 V3 V4 V5 V6 V7 V8 Weight ddb ddth d1st s1 V9 V10 V11 V12 s2 V9 V10 V11 V12
+ * < ncovcol=6 > nqv=2 (V7 V8) dv dv dv qtv dv dv dvv qtv
+ * ntv=3 nqtv=1
+ * cptcovn number of covariates (not including constant and age) = # of + plus 1 = 10+1=11
+ * For time varying covariate, quanti or dummies
+ * cotqvar[wav][iv(1 to nqtv)][i]= [1][12][i]=(V12) quanti
+ * cotvar[wav][ntv+iv][i]= [3+(1 to nqtv)][i]=(V12) quanti
+ * cotvar[wav][iv(1 to ntv)][i]= [1][1][i]=(V9) dummies at wav 1
+ * cotvar[wav][iv(1 to ntv)][i]= [1][2][i]=(V10) dummies at wav 1
+ * covar[k,i], value of kth fixed covariate dummy or quanti :
+ * covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)
+ * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 + V9 + V9*age + V10
+ * k= 1 2 3 4 5 6 7 8 9 10 11
+ */
+/* According to the model, more columns can be added to covar by the product of covariates */
/* ncovcol=1(Males=0 Females=1) nqv=1(raedyrs) ntv=2(withoutiadl=0 withiadl=1, witoutadl=0 withoutadl=1) nqtv=1(bmi) nlstate=3 ndeath=1
# States 1=Coresidence, 2 Living alone, 3 Institution
# V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi
*/
-/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
-/*k 1 2 3 4 5 6 7 8 9 */
-/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */
-/* Tndvar[k] 1 2 3 4 5 */
-/*TDvar 4 3 6 7 1 */ /* For outputs only; combination of dummies fixed or varying */
-/* Tns[k] 1 2 2 4 5 */ /* Number of single cova */
-/* TvarsD[k] 1 2 3 */ /* Number of single dummy cova */
-/* TvarsDind 2 3 9 */ /* position K of single dummy cova */
-/* TvarsQ[k] 1 2 */ /* Number of single quantitative cova */
-/* TvarsQind 1 6 */ /* position K of single quantitative cova */
-/* Tprod[i]=k 4 7 */
-/* Tage[i]=k 5 8 */
-/* */
+/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+/* k 1 2 3 4 5 6 7 8 9 */
+/*Typevar[k]= 0 0 0 2 1 0 2 1 0 *//*0 for simple covariate (dummy, quantitative,*/
+ /* fixed or varying), 1 for age product, 2 for*/
+ /* product */
+/*Dummy[k]= 1 0 0 1 3 1 1 2 0 *//*Dummy[k] 0=dummy (0 1), 1 quantitative */
+ /*(single or product without age), 2 dummy*/
+ /* with age product, 3 quant with age product*/
+/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */
+/* nsd 1 2 3 */ /* Counting single dummies covar fixed or tv */
+/*TvarsD[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/
+/*TvarsDind[k] 2 3 9 */ /* position K of single dummy cova */
+/* nsq 1 2 */ /* Counting single quantit tv */
+/* TvarsQ[k] 5 2 */ /* Number of single quantitative cova */
+/* TvarsQind 1 6 */ /* position K of single quantitative cova */
+/* Tprod[i]=k 1 2 */ /* Position in model of the ith prod without age */
+/* cptcovage 1 2 */ /* Counting cov*age in the model equation */
+/* Tage[cptcovage]=k 5 8 */ /* Position in the model of ith cov*age */
+/* Tvard[1][1]@4={4,3,1,2} V4*V3 V1*V2 */ /* Position in model of the ith prod without age */
+/* TvarF TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */
+/* TvarFind; /**< TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */
/* Type */
/* V 1 2 3 4 5 */
/* F F V V V */
double fp,fptt;
double *xits;
int niterf, itmp;
-#ifdef LINMINORIGINAL
-#else
-
- flatdir=ivector(1,n);
- for (j=1;j<=n;j++) flatdir[j]=0;
-#endif
pt=vector(1,n);
ptt=vector(1,n);
/* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */
/* New value of last point Pn is not computed, P(n-1) */
- for(j=1;j<=n;j++) {
- if(flatdir[j] >0){
- printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- }
- /* printf("\n"); */
- /* fprintf(ficlog,"\n"); */
+ for(j=1;j<=n;j++) {
+ if(flatdir[j] >0){
+ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
}
+ /* printf("\n"); */
+ /* fprintf(ficlog,"\n"); */
+ }
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
}
#endif
-#ifdef LINMINORIGINAL
-#else
- free_ivector(flatdir,1,n);
-#endif
free_vector(xit,1,n);
free_vector(xits,1,n);
free_vector(ptt,1,n);
}
printf("\n");
fprintf(ficlog,"\n");
+#ifdef FLATSUP
+ free_vector(xit,1,n);
+ free_vector(xits,1,n);
+ free_vector(ptt,1,n);
+ free_vector(pt,1,n);
+ return;
+#endif
}
#endif
printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
newm=savm;
/* Covariates have to be included here again */
cov[2]=agefin;
- if(nagesqr==1)
- cov[3]= agefin*agefin;;
+ if(nagesqr==1){
+ cov[3]= agefin*agefin;
+ }
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)];
+ /* cov[++k1]=nbcode[TvarsD[k]][codtabm(ij,k)]; */
/* printf("prevalim 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)); */
}
for (k=1; k<=nsq;k++) { /* For single varying covariates only */
/* Here comes the value of quantitative after renumbering k with single quantitative covariates */
- cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
+ cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
+ /* cov[++k1]=Tqresult[nres][k]; */
/* printf("prevalim Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
}
for (k=1; k<=cptcovage;k++){ /* For product with age */
- if(Dummy[Tvar[Tage[k]]]){
+ if(Dummy[Tage[k]]==2){ /* dummy with age */
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- } else{
- cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
+ /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */
+ } else if(Dummy[Tage[k]]==3){ /* quantitative with age */
+ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
+ /* cov[++k1]=Tqresult[nres][k]; */
}
/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
}
if(Dummy[Tvard[k][1]==0]){
if(Dummy[Tvard[k][2]==0]){
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+ /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
}else{
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
+ /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */
}
}else{
if(Dummy[Tvard[k][2]==0]){
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
+ /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */
}else{
cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]];
+ /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; */
}
}
}
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
/* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
- /* age and covariate values of ij are in 'cov' */
+ /* age and covariate values of ij are in 'cov' */
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */
savm=oldm;
/* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */
/* Covariates have to be included here again */
cov[2]=agefin;
- if(nagesqr==1)
+ if(nagesqr==1){
cov[3]= agefin*agefin;;
+ }
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)];
/* /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
for (k=1; k<=cptcovage;k++){ /* For product with age */
- if(Dummy[Tvar[Tage[k]]]){
+ /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ ERROR ???*/
+ if(Dummy[Tage[k]]== 2){ /* dummy with age */
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- } else{
- cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
+ } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
+ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
}
/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
}
cov[1]=1.;
agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
cov[2]=agexact;
- if(nagesqr==1)
+ if(nagesqr==1){
cov[3]= agexact*agexact;
+ }
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 */
+/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
+ /* codtabm(ij,k) (1 & (ij-1) >> (k-1))+1 */
+/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+/* k 1 2 3 4 5 6 7 8 9 */
+/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */
+/* nsd 1 2 3 */ /* Counting single dummies covar fixed or tv */
+/*TvarsD[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/
+/*TvarsDind[k] 2 3 9 */ /* position K of single dummy cova */
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
/* printf("hpxij 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)); */
}
for (k=1; k<=nsq;k++) { /* For single varying covariates only */
/* Here comes the value of quantitative after renumbering k with single quantitative covariates */
- cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
+ cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
/* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
}
- for (k=1; k<=cptcovage;k++){
- if(Dummy[Tvar[Tage[k]]]){
+ for (k=1; k<=cptcovage;k++){ /* For product with age V1+V1*age +V4 +age*V3 */
+ /* 1+2 Tage[1]=2 TVar[2]=1 Dummy[2]=2, Tage[2]=4 TVar[4]=3 Dummy[4]=3 quant*/
+ /* */
+ if(Dummy[Tage[k]]== 2){ /* dummy with age */
+ /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ */
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- } else{
- cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
+ } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
+ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
}
/* printf("hPxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
}
- for (k=1; k<=cptcovprod;k++){ /* */
+ for (k=1; k<=cptcovprod;k++){ /* For product without age */
/* printf("hPxij Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
- 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,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
+ if(Dummy[Tvard[k][1]==0]){
+ if(Dummy[Tvard[k][2]==0]){
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+ }else{
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
+ }
+ }else{
+ if(Dummy[Tvard[k][2]==0]){
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
+ }else{
+ cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]];
+ }
+ }
}
/* for (k=1; k<=cptcovn;k++) */
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
- /* right multiplication of oldm by the current matrix */
+ /* right multiplication of oldm by the current matrix */
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
pmij(pmmij,cov,ncovmodel,x,nlstate));
/* if((int)age == 70){ */
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
/* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
}
- for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 */
- if(Dummy[Tvar[Tage[k]]]){
+ for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 *//* For product with age */
+ /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age error!!!*\/ */
+ if(Dummy[Tage[k]]== 2){ /* dummy with age */
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- } else{
+ } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
}
/* printf("hBxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
*/
ioffset=2+nagesqr ;
/* Fixed */
- for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */
- cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
+ for (k=1; k<=ncovf;k++){ /* For each fixed covariate dummu or quant or prod */
+ /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */
+ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+ /* TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */
+ /* TvarFind; /**< TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */
+ cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/
+ /* V1*V2 (7) TvarFind[2]=7, TvarFind[3]=9 */
}
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]
- is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]
+ is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6
has been calculated etc */
/* For an individual i, wav[i] gives the number of effective waves */
/* We compute the contribution to Likelihood of each effective transition
meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
*/
for(mi=1; mi<= wav[i]-1; mi++){
- for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/
- /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
+ for(k=1; k <= ncovv ; k++){ /* Varying covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3) Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*/
+ /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? */
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
}
for (ii=1;ii<=nlstate+ndeath;ii++)
} /* end of individual */
} else if(mle==2){
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
+ ioffset=2+nagesqr ;
+ for (k=1; k<=ncovf;k++)
+ cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];
for(mi=1; mi<= wav[i]-1; mi++){
+ for(k=1; k <= ncovv ; k++){
+ cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
+ }
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);
void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))
{
- int i,j, iter=0;
+ int i,j,k, jk, jkk=0, iter=0;
double **xi;
double fret;
double fretone; /* Only one call to likelihood */
if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
fprintf(ficrespow,"\n");
#ifdef POWELL
+#ifdef LINMINORIGINAL
+#else /* LINMINORIGINAL */
+
+ flatdir=ivector(1,npar);
+ for (j=1;j<=npar;j++) flatdir[j]=0;
+#endif /*LINMINORIGINAL */
+
+#ifdef FLATSUP
+ powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
+ /* reorganizing p by suppressing flat directions */
+ for(i=1, jk=1; i <=nlstate; i++){
+ for(k=1; k <=(nlstate+ndeath); k++){
+ if (k != i) {
+ printf("%d%d flatdir[%d]=%d",i,k,jk, flatdir[jk]);
+ if(flatdir[jk]==1){
+ printf(" To be skipped %d%d flatdir[%d]=%d ",i,k,jk, flatdir[jk]);
+ }
+ for(j=1; j <=ncovmodel; j++){
+ printf("%12.7f ",p[jk]);
+ jk++;
+ }
+ printf("\n");
+ }
+ }
+ }
+/* skipping */
+ /* for(i=1, jk=1, jkk=1;(flatdir[jk]==0)&& (i <=nlstate); i++){ */
+ for(i=1, jk=1, jkk=1;i <=nlstate; i++){
+ for(k=1; k <=(nlstate+ndeath); k++){
+ if (k != i) {
+ printf("%d%d flatdir[%d]=%d",i,k,jk, flatdir[jk]);
+ if(flatdir[jk]==1){
+ printf(" To be skipped %d%d flatdir[%d]=%d jk=%d p[%d] ",i,k,jk, flatdir[jk],jk, jk);
+ for(j=1; j <=ncovmodel; jk++,j++){
+ printf(" p[%d]=%12.7f",jk, p[jk]);
+ /*q[jjk]=p[jk];*/
+ }
+ }else{
+ printf(" To be kept %d%d flatdir[%d]=%d jk=%d q[%d]=p[%d] ",i,k,jk, flatdir[jk],jk, jkk, jk);
+ for(j=1; j <=ncovmodel; jk++,jkk++,j++){
+ printf(" p[%d]=%12.7f=q[%d]",jk, p[jk],jkk);
+ /*q[jjk]=p[jk];*/
+ }
+ }
+ printf("\n");
+ }
+ fflush(stdout);
+ }
+ }
+ powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
+#else /* FLATSUP */
powell(p,xi,npar,ftol,&iter,&fret,func);
-#endif
+#endif /* FLATSUP */
+
+#ifdef LINMINORIGINAL
+#else
+ free_ivector(flatdir,1,npar);
+#endif /* LINMINORIGINAL*/
+#endif /* POWELL */
#ifdef NLOPT
#ifdef NEWUOA
}
nlopt_destroy(opt);
#endif
+#ifdef FLATSUP
+ /* npared = npar -flatd/ncovmodel; */
+ /* xired= matrix(1,npared,1,npared); */
+ /* paramred= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */
+ /* powell(pred,xired,npared,ftol,&iter,&fret,flatdir,func); */
+ /* free_matrix(xire,1,npared,1,npared); */
+#else /* FLATSUP */
+#endif /* FLATSUP */
free_matrix(xi,1,npar,1,npar);
fclose(ficrespow);
printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
- fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);
+ fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies (weight=%d) and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm, weightopt);
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
exit(70);
} else{
fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
-<hr size=\"2\" color=\"#EC5E5E\"> \n \
+,<hr size=\"2\" color=\"#EC5E5E\"> \n \
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
- fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>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) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
+ fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>(weight=%d) 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) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr,weightopt);
y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
/* } */
} /* end bool */
} /* end iind = 1 to imx */
- /* prop[s][age] is feeded for any initial and valid live state as well as
+ /* prop[s][age] is fed for any initial and valid live state as well as
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
}
}
-
+ if((int)age ==50){
+ printf(" age=%d cij=%d nres=%d varhe[%d][%d]=%f ",(int)age, cij, nres, 1,2,varhe[1][2]);
+ }
/* Computing expectancies */
hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres);
for(i=1; i<=nlstate;i++)
fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ /* for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); */
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
fprintf(ficresprobcor, "\n#********** Variable ");
*/
/* nbcode[1][1]=0 nbcode[1][2]=1;*/
}
- /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- for (k=1; k<=cptcovage;k++) cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
+ /* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] */
+ /*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ for (k=1; k<=cptcovage;k++)
+ cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
for (k=1; k<=cptcovprod;k++)
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \
double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){
int jj1, k1, i1, cpt, k4, nres;
-
+ /* In fact some results are already printed in fichtm which is open */
fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
<li><a href='#secondorder'>Result files (second order (variance)</a>\n \
</ul>");
- fprintf(fichtm,"<ul><li> model=1+age+%s\n \
-</ul>", model);
+/* fprintf(fichtm,"<ul><li> model=1+age+%s\n \ */
+/* </ul>", model); */
fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
}
if(lval <-1 || lval >1){
printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
- Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+ Should be a value of %d(nth) covariate of wave %d (0 should be the value for the reference and 1\n \
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
For example, for multinomial values like 1, 2 and 3,\n \
build V1=0 V2=0 for the reference value (1),\n \
V1=1 V2=0 for (2) \n \
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
output of IMaCh is often meaningless.\n \
- Exiting.\n",lval,linei, i,line,j);
+ Exiting.\n",lval,linei, i,line,iv,j);
fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
- Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+ Should be a value of %d(nth) covariate of wave %d (0 should be the value for the reference and 1\n \
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
For example, for multinomial values like 1, 2 and 3,\n \
build V1=0 V2=0 for the reference value (1),\n \
V1=1 V2=0 for (2) \n \
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
output of IMaCh is often meaningless.\n \
- Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
+ Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog);
return 1;
}
cotvar[j][iv][i]=(double)(lval);
* - cptcovs number of simple covariates
* - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
* which is a new column after the 9 (ncovcol) variables.
- * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual
+ * - if k is a product Vn*Vm, covar[k][i] is filled with correct values for each individual
* - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
* Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
* - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
*/
+/* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
{
int i, j, k, ks, v;
int j1, k1, k2, k3, k4;
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2
* k= 1 2 3 4 5 6 7 8 9 10 11 12
* Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8
- * p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6}
+ * p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6}
* p Tprod[1]@2={ 6, 5}
*p Tvard[1][1]@4= {7, 8, 5, 6}
* covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
- *How to reorganize?
+ *How to reorganize? Tvars(orted)
* Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age
* Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6}
* {2, 1, 4, 8, 5, 6, 3, 7}
Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;
}
cptcovage=0;
- for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
- cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
- modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */
- if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
+ for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */
+ cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right
+ modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */ /* <model> "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */
+ if (nbocc(modelsav,'+')==0)
+ strcpy(strb,modelsav); /* and analyzes it */
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
/*scanf("%d",i);*/
- if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
- cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
+ if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V5*age+ V4+V3*age strb=V3*age */
+ cutl(strc,strd,strb,'*'); /**< k=1 strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
/* covar is not filled and then is empty */
cptcovprod--;
cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
- Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
+ Tvar[k]=atoi(stre); /* V2+V1+V5*age+V4+V3*age Tvar[5]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
Typevar[k]=1; /* 1 for age product */
- cptcovage++; /* Sums the number of covariates which include age as a product */
- Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
+ cptcovage++; /* Counts the number of covariates which include age as a product */
+ Tage[cptcovage]=k; /* V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
/*printf("stre=%s ", stre);*/
} else if (strcmp(strd,"age")==0) { /* or age*Vn */
cptcovprod--;
Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
because this model-covariate is a construction we invent a new column
which is after existing variables ncovcol+nqv+ntv+nqtv + k1
- If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
- Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
+ If already ncovcol=4 and model=V2 + V1 +V1*V4 +age*V3 +V3*V2
+ thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
+ Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */
Typevar[k]=2; /* 2 for double fixed dummy covariates */
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */
- Tposprod[k]=k1; /* Tpsprod[3]=1, Tposprod[2]=5 */
+ Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */
Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */
}
} /* End age is not in the model */
} /* End if model includes a product */
- else { /* no more sum */
+ else { /* not a product */
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
/* scanf("%d",i);*/
cutl(strd,strc,strb,'V');
model= V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place
k = 1 2 3 4 5 6 7 8 9
Tvar[k]= 5 4 3 1+1+2+1+1=6 5 2 7 1 5
- Typevar[k]= 0 0 0 2 1 0 2 1 1
+ Typevar[k]= 0 0 0 2 1 0 2 1 0
Fixed[k] 1 1 1 1 3 0 0 or 2 2 3
Dummy[k] 1 0 0 0 3 1 1 2 3
Tmodelind[combination of covar]=k;
modell[k].subtype= VQ;
ncovv++; /* Only simple time varying variables */
nsq++;
- TvarsQ[nsq]=Tvar[k];
+ TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */
TvarsQind[nsq]=k;
TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
double dum=0.; /* Dummy variable */
double ***p3mat;
/* double ***mobaverage; */
+ double wald;
char line[MAXLINE];
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
- printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+ printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); /* Printing model equation */
fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+
+ printf("#model= 1 + age ");
+ fprintf(ficres,"#model= 1 + age ");
+ fprintf(ficlog,"#model= 1 + age ");
+ fprintf(fichtm,"\n<ul><li> model=1+age+%s\n \
+</ul>", model);
+
+ fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">\n");
+ fprintf(fichtm, "<tr><th>Model=</th><th>1</th><th>+ age</th>");
+ if(nagesqr==1){
+ printf(" + age*age ");
+ fprintf(ficres," + age*age ");
+ fprintf(ficlog," + age*age ");
+ fprintf(fichtm, "<th>+ age*age</th>");
+ }
+ for(j=1;j <=ncovmodel-2;j++){
+ if(Typevar[j]==0) {
+ printf(" + V%d ",Tvar[j]);
+ fprintf(ficres," + V%d ",Tvar[j]);
+ fprintf(ficlog," + V%d ",Tvar[j]);
+ fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
+ }else if(Typevar[j]==1) {
+ printf(" + V%d*age ",Tvar[j]);
+ fprintf(ficres," + V%d*age ",Tvar[j]);
+ fprintf(ficlog," + V%d*age ",Tvar[j]);
+ fprintf(fichtm, "<th>+ V%d*age</th>",Tvar[j]);
+ }else if(Typevar[j]==2) {
+ printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+ fprintf(ficres," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+ fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+ fprintf(fichtm, "<th>+ V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+ }
+ }
+ printf("\n");
+ fprintf(ficres,"\n");
+ fprintf(ficlog,"\n");
+ fprintf(fichtm, "</tr>");
+ fprintf(fichtm, "\n");
+
+
for(i=1,jk=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
+ fprintf(fichtm, "<tr>");
printf("%d%d ",i,k);
fprintf(ficlog,"%d%d ",i,k);
fprintf(ficres,"%1d%1d ",i,k);
+ fprintf(fichtm, "<td>%1d%1d</td>",i,k);
for(j=1; j <=ncovmodel; j++){
printf("%12.7f ",p[jk]);
fprintf(ficlog,"%12.7f ",p[jk]);
fprintf(ficres,"%12.7f ",p[jk]);
+ fprintf(fichtm, "<td>%12.7f</td>",p[jk]);
jk++;
}
printf("\n");
fprintf(ficlog,"\n");
fprintf(ficres,"\n");
+ fprintf(fichtm, "</tr>\n");
}
}
}
+ /* fprintf(fichtm,"</tr>\n"); */
+ fprintf(fichtm,"</table>\n");
+ fprintf(fichtm, "\n");
+
if(mle != 0){
/* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */
ftolhess=ftol; /* Usually correct */
hesscov(matcov, hess, p, npar, delti, ftolhess, func);
printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
+ fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">");
+ fprintf(fichtm, "\n<tr><th>Model=</th><th>1</th><th>+ age</th>");
+ if(nagesqr==1){
+ printf(" + age*age ");
+ fprintf(ficres," + age*age ");
+ fprintf(ficlog," + age*age ");
+ fprintf(fichtm, "<th>+ age*age</th>");
+ }
+ for(j=1;j <=ncovmodel-2;j++){
+ if(Typevar[j]==0) {
+ printf(" + V%d ",Tvar[j]);
+ fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
+ }else if(Typevar[j]==1) {
+ printf(" + V%d*age ",Tvar[j]);
+ fprintf(fichtm, "<th>+ V%d*age</th>",Tvar[j]);
+ }else if(Typevar[j]==2) {
+ fprintf(fichtm, "<th>+ V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+ }
+ }
+ fprintf(fichtm, "</tr>\n");
+
for(i=1,jk=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
+ fprintf(fichtm, "<tr valign=top>");
printf("%d%d ",i,k);
fprintf(ficlog,"%d%d ",i,k);
+ fprintf(fichtm, "<td>%1d%1d</td>",i,k);
for(j=1; j <=ncovmodel; j++){
- printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
- fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+ wald=p[jk]/sqrt(matcov[jk][jk]);
+ printf("%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+ fprintf(ficlog,"%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+ if(fabs(wald) > 1.96){
+ fprintf(fichtm, "<td><b>%12.7f</b> (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));
+ fprintf(fichtm,"<b>W=%8.3f</b></br>",wald);
+ }else{
+ fprintf(fichtm, "<td>%12.7f (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));
+ fprintf(fichtm,"W=%8.3f</br>",wald);
+ }
+ fprintf(fichtm,"[%12.7f;%12.7f]</br></td>", p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
jk++;
}
printf("\n");
fprintf(ficlog,"\n");
+ fprintf(fichtm, "</tr>\n");
}
}
}
} /* end of hesscov and Wald tests */
+ fprintf(fichtm,"</table>\n");
/* */
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
git://henry.ined.fr / .git / commitdiff