/* $Id$
$State$
$Log$
+ Revision 1.233 2016/08/23 07:40:50 brouard
+ Summary: not working
+
Revision 1.232 2016/08/22 14:20:21 brouard
Summary: not working
int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */
int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */
int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */
-
+int nsd=0; /**< Total number of single dummy variables (output) */
+int nsq=0; /**< Total number of single quantitative variables (output) */
int ncoveff=0; /* Total number of effective fixed dummy covariates in the model */
int nqfveff=0; /**< nqfveff Number of Quantitative Fixed Variables Effective */
int ntveff=0; /**< ntveff number of effective time varying variables */
FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];
char title[MAXLINE];
+char model[MAXLINE]; /**< The model line */
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH];
char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH];
double ***cotqvar; /* Time varying quantitative covariate itqv */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+/* 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 */
+/* */
+/* Type */
+/* V 1 2 3 4 5 */
+/* F F V V V */
+/* D Q D D Q */
+/* */
+int *TvarsD;
+int *TvarsDind;
+int *TvarsQ;
+int *TvarsQind;
+
+/* int *TDvar; /\**< TDvar[1]=4, TDvarF[2]=3, TDvar[3]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */
int *TvarF; /**< TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
int *TvarFind; /**< TvarFind[1]=6, TvarFind[2]=7, Tvarind[3]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
int *TvarV; /**< TvarV[1]=Tvar[1]=5, TvarV[2]=Tvar[2]=4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
i=0;
lg=strlen(s);
for(i=0; i<= lg; i++) {
- if (s[i] == occ ) j++;
+ if (s[i] == occ ) j++;
}
return j;
}
if (directest < 0.0) { /* Then we use it for new direction */
#endif
#ifdef DEBUGLINMIN
- printf("Before linmin in direction P%d-P0\n",n);
- for (j=1;j<=n;j++) {
- printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0){
- printf("\n");
- fprintf(ficlog,"\n");
- }
- }
+ printf("Before linmin in direction P%d-P0\n",n);
+ for (j=1;j<=n;j++) {
+ printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
+ printf("\n");
+ fprintf(ficlog,"\n");
+ }
+ }
#endif
#ifdef LINMINORIGINAL
- linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+ linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
#else
- linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
- flatdir[i]=flat; /* Function is vanishing in that direction i */
+ linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+ flatdir[i]=flat; /* Function is vanishing in that direction i */
#endif
-
+
#ifdef DEBUGLINMIN
- for (j=1;j<=n;j++) {
- printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0){
- printf("\n");
- fprintf(ficlog,"\n");
- }
- }
+ for (j=1;j<=n;j++) {
+ printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
+ printf("\n");
+ fprintf(ficlog,"\n");
+ }
+ }
#endif
- for (j=1;j<=n;j++) {
- xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
- xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
- }
+ for (j=1;j<=n;j++) {
+ xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
+ xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
+ }
#ifdef LINMINORIGINAL
#else
- for (j=1, flatd=0;j<=n;j++) {
- if(flatdir[j]>0)
- flatd++;
- }
- if(flatd >0){
- printf("%d flat directions\n",flatd);
- fprintf(ficlog,"%d flat directions\n",flatd);
- for (j=1;j<=n;j++) {
- if(flatdir[j]>0){
- printf("%d ",j);
- fprintf(ficlog,"%d ",j);
- }
- }
- printf("\n");
- fprintf(ficlog,"\n");
- }
+ for (j=1, flatd=0;j<=n;j++) {
+ if(flatdir[j]>0)
+ flatd++;
+ }
+ if(flatd >0){
+ printf("%d flat directions\n",flatd);
+ fprintf(ficlog,"%d flat directions\n",flatd);
+ for (j=1;j<=n;j++) {
+ if(flatdir[j]>0){
+ printf("%d ",j);
+ fprintf(ficlog,"%d ",j);
+ }
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
+ }
#endif
- printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
- fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
-
+ printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+ fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+
#ifdef DEBUG
- printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
- fprintf(ficlog,"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(" %lf",xit[j]);
- fprintf(ficlog," %lf",xit[j]);
- }
- printf("\n");
- fprintf(ficlog,"\n");
+ printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+ fprintf(ficlog,"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(" %lf",xit[j]);
+ fprintf(ficlog," %lf",xit[j]);
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
#endif
} /* end of t or directest negative */
#ifdef POWELLNOF3INFF1TEST
#else
- } /* end if (fptt < fp) */
+ } /* end if (fptt < fp) */
#endif
#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */
- } /*NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done */
+ } /*NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done */
#else
#endif
- } /* loop iteration */
+ } /* loop iteration */
}
-
+
/**** Prevalence limit (stable or period prevalence) ****************/
-
+
double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
-{
- /* Computes the prevalence limit in each live state at age x and for covariate ij by left multiplying the unit
- matrix by transitions matrix until convergence is reached with precision ftolpl */
+ {
+ /* Computes the prevalence limit in each live state at age x and for covariate combiation ij by left multiplying the unit
+ matrix by transitions matrix until convergence is reached with precision ftolpl */
/* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */
/* Wx is row vector: population in state 1, population in state 2, population dead */
/* or prevalence in state 1, prevalence in state 2, 0 */
/* {0.51571254859325999, 0.4842874514067399, */
/* 0.51326036147820708, 0.48673963852179264} */
/* If we start from prlim again, prlim tends to a constant matrix */
-
+
int i, ii,j,k;
double *min, *max, *meandiff, maxmax,sumnew=0.;
/* double **matprod2(); */ /* test */
cov[2]=agefin;
if(nagesqr==1)
cov[3]= agefin*agefin;;
- for (k=1; k<=cptcovn;k++) {
- /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
- /* Here comes the value of the covariate 'ij' */
- cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
- /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
+ for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
+ /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
+ cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
+ printf("prevalim ij=%d k=%d TvarsD[%d]=%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]]=qselvar[k]; */
+ printf("prevalim ij=%d k=%d TvarsQind[%d]=%d \n",ij,k,k,TvarsQind[k]);
}
/*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
/* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
- for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
- for (k=1; k<=cptcovprod;k++) /* Useless */
- /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+ for (k=1; k<=cptcovage;k++){
+ if(Dummy[Tvar[Tage[k]]]){
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ } else{
+ ;
+ /* cov[2+nagesqr+Tage[k]]=qselvar[k]; */
+ }
+ printf("prevalim Age ij=%d k=%d Tage[%d]=%d \n",ij,k,k,Tage[k]);
+ }
+ for (k=1; k<=cptcovprod;k++){ /* */
+ printf("prevalim Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=%d, Tvard[%d][2]=%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)];
-
+ }
/*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]);*/
}
for(j=1; j<=nlstate; j++){
for(i=1;i<=nlstate;i++){
- /* bprlim[i][j]= newm[i][j]/(1-sumnew); */
- bprlim[i][j]= newm[i][j];
- max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */
- min[i]=FMIN(min[i],bprlim[i][j]);
+ /* bprlim[i][j]= newm[i][j]/(1-sumnew); */
+ bprlim[i][j]= newm[i][j];
+ max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */
+ min[i]=FMIN(min[i],bprlim[i][j]);
}
}
/*double t34;*/
int i,j, nc, ii, jj;
- for(i=1; i<= nlstate; i++){
- for(j=1; j<i;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += param[i][j][nc]*cov[nc];*/
- lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
- /* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
- }
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- /* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
- }
- for(j=i+1; j<=nlstate+ndeath;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
- lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
- /* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
- }
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- }
- }
-
- for(i=1; i<= nlstate; i++){
- s1=0;
- for(j=1; j<i; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- for(j=i+1; j<=nlstate+ndeath; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
- ps[i][i]=1./(s1+1.);
- /* Computing other pijs */
- for(j=1; j<i; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- for(j=i+1; j<=nlstate+ndeath; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
- } /* end i */
-
- for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
- for(jj=1; jj<= nlstate+ndeath; jj++){
- ps[ii][jj]=0;
- ps[ii][ii]=1;
- }
- }
- /* Added for backcast */ /* Transposed matrix too */
- for(jj=1; jj<= nlstate+ndeath; jj++){
- s1=0.;
- for(ii=1; ii<= nlstate+ndeath; ii++){
- s1+=ps[ii][jj];
- }
- for(ii=1; ii<= nlstate; ii++){
- ps[ii][jj]=ps[ii][jj]/s1;
- }
- }
- /* Transposition */
- for(jj=1; jj<= nlstate+ndeath; jj++){
- for(ii=jj; ii<= nlstate+ndeath; ii++){
- s1=ps[ii][jj];
- ps[ii][jj]=ps[jj][ii];
- ps[jj][ii]=s1;
- }
- }
- /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
- /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
- /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
- /* } */
- /* printf("\n "); */
- /* } */
- /* printf("\n ");printf("%lf ",cov[2]);*/
- /*
- for(i=1; i<= npar; i++) printf("%f ",x[i]);
- goto end;*/
- return ps;
+ for(i=1; i<= nlstate; i++){
+ for(j=1; j<i;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += param[i][j][nc]*cov[nc];*/
+ lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+ /* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ /* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ for(j=i+1; j<=nlstate+ndeath;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+ lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+ /* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ }
+ }
+
+ for(i=1; i<= nlstate; i++){
+ s1=0;
+ for(j=1; j<i; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ for(j=i+1; j<=nlstate+ndeath; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+ ps[i][i]=1./(s1+1.);
+ /* Computing other pijs */
+ for(j=1; j<i; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ for(j=i+1; j<=nlstate+ndeath; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+ } /* end i */
+
+ for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ ps[ii][jj]=0;
+ ps[ii][ii]=1;
+ }
+ }
+ /* Added for backcast */ /* Transposed matrix too */
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ s1=0.;
+ for(ii=1; ii<= nlstate+ndeath; ii++){
+ s1+=ps[ii][jj];
+ }
+ for(ii=1; ii<= nlstate; ii++){
+ ps[ii][jj]=ps[ii][jj]/s1;
+ }
+ }
+ /* Transposition */
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ for(ii=jj; ii<= nlstate+ndeath; ii++){
+ s1=ps[ii][jj];
+ ps[ii][jj]=ps[jj][ii];
+ ps[jj][ii]=s1;
+ }
+ }
+ /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+ /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
+ /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+ /* } */
+ /* printf("\n "); */
+ /* } */
+ /* printf("\n ");printf("%lf ",cov[2]);*/
+ /*
+ for(i=1; i<= npar; i++) printf("%f ",x[i]);
+ goto end;*/
+ return ps;
}
*/
ioffset=2+nagesqr+cptcovage;
/* 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++){ /* 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)*/
+ }
/* 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]
has been calculated etc */
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 (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
- oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- savm[ii][j]=(ii==j ? 1.0 : 0.0);
- }
- for(d=0; d<dh[mi][i]; d++){
- newm=savm;
- agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
- cov[2]=agexact;
- if(nagesqr==1)
- cov[3]= agexact*agexact; /* Should be changed here */
- for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
- }
- out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- savm=oldm;
- oldm=newm;
- } /* end mult */
-
- /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
- /* But now since version 0.9 we anticipate for bias at large stepm.
- * If stepm is larger than one month (smallest stepm) and if the exact delay
- * (in months) between two waves is not a multiple of stepm, we rounded to
- * the nearest (and in case of equal distance, to the lowest) interval but now
- * we keep into memory the bias bh[mi][i] and also the previous matrix product
- * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
- * probability in order to take into account the bias as a fraction of the way
+ 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 (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+ for(d=0; d<dh[mi][i]; d++){
+ newm=savm;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact; /* Should be changed here */
+ for (kk=1; kk<=cptcovage;kk++) {
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+ }
+ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+ 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+ savm=oldm;
+ oldm=newm;
+ } /* end mult */
+
+ /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+ /* But now since version 0.9 we anticipate for bias at large stepm.
+ * If stepm is larger than one month (smallest stepm) and if the exact delay
+ * (in months) between two waves is not a multiple of stepm, we rounded to
+ * the nearest (and in case of equal distance, to the lowest) interval but now
+ * we keep into memory the bias bh[mi][i] and also the previous matrix product
+ * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
+ * probability in order to take into account the bias as a fraction of the way
* from savm to out if bh is negative or even beyond if bh is positive. bh varies
* -stepm/2 to stepm/2 .
* For stepm=1 the results are the same as for previous versions of Imach.
* For stepm > 1 the results are less biased than in previous versions.
*/
- s1=s[mw[mi][i]][i];
- s2=s[mw[mi+1][i]][i];
- bbh=(double)bh[mi][i]/(double)stepm;
- /* bias bh is positive if real duration
- * is higher than the multiple of stepm and negative otherwise.
- */
- /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
- if( s2 > nlstate){
- /* i.e. if s2 is a death state and if the date of death is known
- then the contribution to the likelihood is the probability to
- die between last step unit time and current step unit time,
- which is also equal to probability to die before dh
- minus probability to die before dh-stepm .
- In version up to 0.92 likelihood was computed
- as if date of death was unknown. Death was treated as any other
- health state: the date of the interview describes the actual state
- and not the date of a change in health state. The former idea was
- to consider that at each interview the state was recorded
- (healthy, disable or death) and IMaCh was corrected; but when we
- introduced the exact date of death then we should have modified
- the contribution of an exact death to the likelihood. This new
- contribution is smaller and very dependent of the step unit
- stepm. It is no more the probability to die between last interview
- and month of death but the probability to survive from last
- interview up to one month before death multiplied by the
- probability to die within a month. Thanks to Chris
- Jackson for correcting this bug. Former versions increased
- mortality artificially. The bad side is that we add another loop
- which slows down the processing. The difference can be up to 10%
- lower mortality.
- */
- /* If, at the beginning of the maximization mostly, the
- cumulative probability or probability to be dead is
- constant (ie = 1) over time d, the difference is equal to
- 0. out[s1][3] = savm[s1][3]: probability, being at state
- s1 at precedent wave, to be dead a month before current
- wave is equal to probability, being at state s1 at
- precedent wave, to be dead at mont of the current
- wave. Then the observed probability (that this person died)
- is null according to current estimated parameter. In fact,
- it should be very low but not zero otherwise the log go to
- infinity.
- */
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ /* bias bh is positive if real duration
+ * is higher than the multiple of stepm and negative otherwise.
+ */
+ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
+ if( s2 > nlstate){
+ /* i.e. if s2 is a death state and if the date of death is known
+ then the contribution to the likelihood is the probability to
+ die between last step unit time and current step unit time,
+ which is also equal to probability to die before dh
+ minus probability to die before dh-stepm .
+ In version up to 0.92 likelihood was computed
+ as if date of death was unknown. Death was treated as any other
+ health state: the date of the interview describes the actual state
+ and not the date of a change in health state. The former idea was
+ to consider that at each interview the state was recorded
+ (healthy, disable or death) and IMaCh was corrected; but when we
+ introduced the exact date of death then we should have modified
+ the contribution of an exact death to the likelihood. This new
+ contribution is smaller and very dependent of the step unit
+ stepm. It is no more the probability to die between last interview
+ and month of death but the probability to survive from last
+ interview up to one month before death multiplied by the
+ probability to die within a month. Thanks to Chris
+ Jackson for correcting this bug. Former versions increased
+ mortality artificially. The bad side is that we add another loop
+ which slows down the processing. The difference can be up to 10%
+ lower mortality.
+ */
+ /* If, at the beginning of the maximization mostly, the
+ cumulative probability or probability to be dead is
+ constant (ie = 1) over time d, the difference is equal to
+ 0. out[s1][3] = savm[s1][3]: probability, being at state
+ s1 at precedent wave, to be dead a month before current
+ wave is equal to probability, being at state s1 at
+ precedent wave, to be dead at mont of the current
+ wave. Then the observed probability (that this person died)
+ is null according to current estimated parameter. In fact,
+ it should be very low but not zero otherwise the log go to
+ infinity.
+ */
/* #ifdef INFINITYORIGINAL */
/* lli=log(out[s1][s2] - savm[s1][s2]); */
/* #else */
for (iind=1; iind<=imx; iind++) { /* For each individual iind */
bool=1;
if(anyvaryingduminmodel==0){ /* If All fixed covariates */
- if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+ if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
/* for (z1=1; z1<= nqfveff; z1++) { */
/* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */
/* } */
- for (z1=1; z1<=cptcoveff; z1++) {
- /* if(Tvaraff[z1] ==-20){ */
- /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
- /* }else if(Tvaraff[z1] ==-10){ */
- /* /\* sumnew+=coqvar[z1][iind]; *\/ */
- /* }else */
- if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
- /* Tests if this individual iind responded to j1 (V4=1 V3=0) */
- bool=0;
- /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",
- bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
- j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
- /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
- } /* Onlyf fixed */
- } /* end z1 */
- } /* cptcovn > 0 */
+ for (z1=1; z1<=cptcoveff; z1++) {
+ /* if(Tvaraff[z1] ==-20){ */
+ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
+ /* }else if(Tvaraff[z1] ==-10){ */
+ /* /\* sumnew+=coqvar[z1][iind]; *\/ */
+ /* }else */
+ if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
+ /* Tests if this individual iind responded to j1 (V4=1 V3=0) */
+ bool=0;
+ /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",
+ bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
+ j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+ /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
+ } /* Onlyf fixed */
+ } /* end z1 */
+ } /* cptcovn > 0 */
} /* end any */
if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */
- /* for(m=firstpass; m<=lastpass; m++){ */
- for(mi=1; mi<wav[iind];mi++){ /* For that wave */
- m=mw[mi][iind];
- if(anyvaryingduminmodel==1){ /* Some are varying covariates */
- for (z1=1; z1<=cptcoveff; z1++) {
- if( Fixed[Tmodelind[z1]]==1){
- iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
- if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */
- bool=0;
- }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
- if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
- bool=0;
- }
- }
- }
- }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */
- /* bool =0 we keep that guy which corresponds to the combination of dummy values */
- if(bool==1){
- /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
- and mw[mi+1][iind]. dh depends on stepm. */
- agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
- ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
- if(m >=firstpass && m <=lastpass){
- k2=anint[m][iind]+(mint[m][iind]/12.);
- /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
- if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */
- if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */
- if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */
- prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
- if (m<lastpass) {
- /* if(s[m][iind]==4 && s[m+1][iind]==4) */
- /* printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
- if(s[m][iind]==-1)
- printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
- freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
- /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
- freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
- }
- } /* end if between passes */
- if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
- dateintsum=dateintsum+k2;
- k2cpt++;
- /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
- }
- } /* end bool 2 */
- } /* end m */
+ /* for(m=firstpass; m<=lastpass; m++){ */
+ for(mi=1; mi<wav[iind];mi++){ /* For that wave */
+ m=mw[mi][iind];
+ if(anyvaryingduminmodel==1){ /* Some are varying covariates */
+ for (z1=1; z1<=cptcoveff; z1++) {
+ if( Fixed[Tmodelind[z1]]==1){
+ iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+ if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */
+ bool=0;
+ }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
+ if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
+ bool=0;
+ }
+ }
+ }
+ }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */
+ /* bool =0 we keep that guy which corresponds to the combination of dummy values */
+ if(bool==1){
+ /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
+ and mw[mi+1][iind]. dh depends on stepm. */
+ agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
+ ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
+ if(m >=firstpass && m <=lastpass){
+ k2=anint[m][iind]+(mint[m][iind]/12.);
+ /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
+ if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */
+ if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */
+ if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */
+ prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
+ if (m<lastpass) {
+ /* if(s[m][iind]==4 && s[m+1][iind]==4) */
+ /* printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
+ if(s[m][iind]==-1)
+ printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
+ freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
+ /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
+ freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
+ }
+ } /* end if between passes */
+ if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
+ dateintsum=dateintsum+k2;
+ k2cpt++;
+ /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
+ }
+ } /* end bool 2 */
+ } /* end m */
} /* end bool */
} /* end iind = 1 to imx */
/* prop[s][age] is feeded for any initial and valid live state as well as
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
for (z1=1; z1<=cptcoveff; z1++){
- fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
}
fprintf(ficresp, "**********\n#");
fprintf(ficresphtm, "**********</h3>\n");
fprintf(ficresphtmfr,"<th>Age</th> ");
for(jk=-1; jk <=nlstate+ndeath; jk++){
for(m=-1; m <=nlstate+ndeath; m++){
- if(jk!=0 && m!=0)
- fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
+ if(jk!=0 && m!=0)
+ fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
}
}
fprintf(ficresphtmfr, "\n");
return (1);
}
-void removespace(char **stri){/*, char stro[]) {*/
+void removefirstspace(char **stri){/*, char stro[]) {*/
char *p1 = *stri, *p2 = *stri;
- do
- while (*p2 == ' ')
- p2++;
- while (*p1++ == *p2++);
- *stri=p1;
+ if (*p2 == ' ')
+ p2++;
+ /* while ((*p1++ = *p2++) !=0) */
+ /* ; */
+ /* do */
+ /* while (*p2 == ' ') */
+ /* p2++; */
+ /* while (*p1++ == *p2++); */
+ *stri=p2;
}
int decoderesult ( char resultline[])
/**< This routine decode one result line and returns the combination # of dummy covariates only **/
{
- int j=0, k=0;
+ int j=0, k=0, k1=0, k2=0, match=0;
char resultsav[MAXLINE];
+ int resultmodel[MAXLINE];
+ int modelresult[MAXLINE];
char stra[80], strb[80], strc[80], strd[80],stre[80];
- removespace(&resultline);
+ removefirstspace(&resultline);
printf("decoderesult:%s\n",resultline);
if (strstr(resultline,"v") !=0){
if (strlen(resultsav) >1){
j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */
}
-
- for(k=1; k<=j;k++){ /* Loop on total covariates of the model */
- cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' '
- resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */
- cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */
+ if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
+ printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
+ fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
+ }
+ for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */
+ if(nbocc(resultsav,'=') >1){
+ cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' '
+ resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */
+ cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */
+ }else
+ cutl(strc,strd,resultsav,'=');
Tvalsel[k]=atof(strc); /* 1 */
-
+
cutl(strc,stre,strd,'V'); /* strd='V4' strc=4 stre='V' */;
Tvarsel[k]=atoi(strc);
/* Typevarsel[k]=1; /\* 1 for age product *\/ */
if (nbocc(stra,'=') >0)
strcpy(resultsav,stra); /* and analyzes it */
}
+ /* Checking if no missing or useless values in comparison of current model needs */
+ for(k1=1; k1<= cptcovt ;k1++){ /* model line */
+ if(Typevar[k1]==0){
+ match=0;
+ for(k2=1; k2 <=j;k2++){
+ if(Tvar[k1]==Tvarsel[k2]) {
+ modelresult[k2]=k1;
+ match=1;
+ break;
+ }
+ }
+ if(match == 0){
+ printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model);
+ }
+ }
+ }
+
+ for(k2=1; k2 <=j;k2++){ /* result line */
+ match=0;
+ for(k1=1; k1<= cptcovt ;k1++){ /* model line */
+ if(Typevar[k1]==0){
+ if(Tvar[k1]==Tvarsel[k2]) {
+ resultmodel[k1]=k2;
+ ++match;
+ }
+ }
+ }
+ if(match == 0){
+ printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model);
+ }else if(match > 1){
+ printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model);
+ }
+ }
+
+ /* We need to deduce which combination number is chosen and save quantitative values */
+
return (0);
}
int selected( int kvar){ /* Selected combination of covariates */
if ((strpt=strstr(model,"age*age")) !=0){
printf(" strpt=%s, model=%s\n",strpt, model);
if(strpt != model){
- printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model);
- fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model); fflush(ficlog);
- return 1;
+ return 1;
}
nagesqr=1;
if (strstr(model,"+age*age") !=0)
- substrchaine(modelsav, model, "+age*age");
+ substrchaine(modelsav, model, "+age*age");
else if (strstr(model,"age*age+") !=0)
- substrchaine(modelsav, model, "age*age+");
+ substrchaine(modelsav, model, "age*age+");
else
- substrchaine(modelsav, model, "age*age");
+ substrchaine(modelsav, model, "age*age");
}else
nagesqr=0;
if (strlen(modelsav) >1){
}
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 '+'
+ 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 */
- /* 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 (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 */
- 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 */
- /*printf("stre=%s ", stre);*/
- } else if (strcmp(strd,"age")==0) { /* or age*Vn */
- cptcovprod--;
- cutl(stre,strb,strc,'V');
- Tvar[k]=atoi(stre);
- Typevar[k]=1; /* 1 for age product */
- cptcovage++;
- Tage[cptcovage]=k;
- } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/
- /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
- cptcovn++;
- cptcovprodnoage++;k1++;
- cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
- 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 */
- 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 */
- 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 */
- /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */
- /* Tvar[cptcovt+k2+1]=Tvard[k1][2]; /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */
+ 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 (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 */
+ 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 */
+ /*printf("stre=%s ", stre);*/
+ } else if (strcmp(strd,"age")==0) { /* or age*Vn */
+ cptcovprod--;
+ cutl(stre,strb,strc,'V');
+ Tvar[k]=atoi(stre);
+ Typevar[k]=1; /* 1 for age product */
+ cptcovage++;
+ Tage[cptcovage]=k;
+ } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/
+ /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+ cptcovn++;
+ cptcovprodnoage++;k1++;
+ cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+ 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 */
+ 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 */
+ 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 */
+ /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */
+ /* Tvar[cptcovt+k2+1]=Tvard[k1][2]; /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */
/*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */
- /* 1 2 3 4 5 | Tvar[5+1)=1, Tvar[7]=2 */
- for (i=1; i<=lastobs;i++){
- /* Computes the new covariate which is a product of
- covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
- covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
- }
- } /* End age is not in the model */
- } /* End if model includes a product */
- else { /* no more sum */
- /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
- /* scanf("%d",i);*/
- cutl(strd,strc,strb,'V');
- ks++; /**< Number of simple covariates dummy or quantitative, fixe or varying */
- cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */
- Tvar[k]=atoi(strd);
- Typevar[k]=0; /* 0 for simple covariates */
- }
- strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
+ /* 1 2 3 4 5 | Tvar[5+1)=1, Tvar[7]=2 */
+ for (i=1; i<=lastobs;i++){
+ /* Computes the new covariate which is a product of
+ covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+ covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+ }
+ } /* End age is not in the model */
+ } /* End if model includes a product */
+ else { /* no more sum */
+ /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
+ /* scanf("%d",i);*/
+ cutl(strd,strc,strb,'V');
+ ks++; /**< Number of simple covariates dummy or quantitative, fixe or varying */
+ cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */
+ Tvar[k]=atoi(strd);
+ Typevar[k]=0; /* 0 for simple covariates */
+ }
+ strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
scanf("%d",i);*/
} /* end of loop + on total covariates */
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
- for(k=1, ncovf=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
- if (Tvar[k] <=ncovcol && (Typevar[k]==0 || Typevar[k]==2)){ /* Simple or product fixed dummy (<=ncovcol) covariates */
+ for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
+ if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
Fixed[k]= 0;
Dummy[k]= 0;
ncoveff++;
ncovf++;
- modell[k].maintype= FTYPE;
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
+ nsd++;
+ modell[k].maintype= FTYPE;
+ TvarsD[nsd]=Tvar[k];
+ TvarsDind[nsd]=k;
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
+ TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+ TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+ }else if( Tvar[k] <=ncovcol && Typevar[k]==2){ /* Product of fixed dummy (<=ncovcol) covariates */
+ Fixed[k]= 0;
+ Dummy[k]= 0;
+ ncoveff++;
+ ncovf++;
+ modell[k].maintype= FTYPE;
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */
Fixed[k]= 0;
Dummy[k]= 1;
nqfveff++;
- modell[k].maintype= FTYPE;
- modell[k].subtype= FQ;
+ modell[k].maintype= FTYPE;
+ modell[k].subtype= FQ;
+ nsq++;
+ TvarsQ[nsq]=Tvar[k];
+ TvarsQind[nsq]=k;
ncovf++;
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
- }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){
+ }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying variables */
Fixed[k]= 1;
Dummy[k]= 0;
ntveff++; /* Only simple time varying dummy variable */
- modell[k].maintype= VTYPE;
- modell[k].subtype= VD;
- ncovv++; /* Only simple time varying variables */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VD;
+ nsd++;
+ TvarsD[nsd]=Tvar[k];
+ TvarsDind[nsd]=k;
+ ncovv++; /* Only simple time varying variables */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */
TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */
printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);
printf("Quasi TmodelInvind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv);
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){ /* Only simple time varying quantitative variable V5*/
- Fixed[k]= 1;
- Dummy[k]= 1;
- nqtveff++;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VQ;
- ncovv++; /* Only simple time varying variables */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ nqtveff++;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VQ;
+ ncovv++; /* Only simple time varying variables */
+ nsq++;
+ TvarsQ[nsq]=Tvar[k];
+ TvarsQind[nsq]=k;
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
- TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
- /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */
- printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv);
+ TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
+ /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */
+ printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv);
printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv);
}else if (Typevar[k] == 1) { /* product with age */
- ncova++;
- TvarA[ncova]=Tvar[k];
- TvarAind[ncova]=k;
+ ncova++;
+ TvarA[ncova]=Tvar[k];
+ TvarAind[ncova]=k;
if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
- Fixed[k]= 2;
- Dummy[k]= 2;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APFD;
- /* ncoveff++; */
+ Fixed[k]= 2;
+ Dummy[k]= 2;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APFD;
+ /* ncoveff++; */
}else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
- Fixed[k]= 2;
- Dummy[k]= 3;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APFQ; /* Product age * fixed quantitative */
- /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */
+ Fixed[k]= 2;
+ Dummy[k]= 3;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APFQ; /* Product age * fixed quantitative */
+ /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */
}else if( Tvar[k] <=ncovcol+nqv+ntv ){
- Fixed[k]= 3;
- Dummy[k]= 2;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APVD; /* Product age * varying dummy */
- /* ntveff++; /\* Only simple time varying dummy variable *\/ */
+ Fixed[k]= 3;
+ Dummy[k]= 2;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APVD; /* Product age * varying dummy */
+ /* ntveff++; /\* Only simple time varying dummy variable *\/ */
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 3;
- Dummy[k]= 3;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APVQ; /* Product age * varying quantitative */
- /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
+ Fixed[k]= 3;
+ Dummy[k]= 3;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APVQ; /* Product age * varying quantitative */
+ /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
}
}else if (Typevar[k] == 2) { /* product without age */
k1=Tposprod[k];
if(Tvard[k1][1] <=ncovcol){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 1;
- Dummy[k]= 0;
- modell[k].maintype= FTYPE;
- modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */
- ncovf++; /* Fixed variables without age */
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv){
- Fixed[k]= 0; /* or 2 ?*/
- Dummy[k]= 1;
- modell[k].maintype= FTYPE;
- modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */
- ncovf++; /* Varying variables without age */
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 0;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ modell[k].maintype= FTYPE;
+ modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */
+ ncovf++; /* Fixed variables without age */
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 0; /* or 2 ?*/
+ Dummy[k]= 1;
+ modell[k].maintype= FTYPE;
+ modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */
+ ncovf++; /* Varying variables without age */
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else if(Tvard[k1][1] <=ncovcol+nqv){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 0; /* or 2 ?*/
- Dummy[k]= 1;
- modell[k].maintype= FTYPE;
- modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */
- ncovf++; /* Fixed variables without age */
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 0; /* or 2 ?*/
+ Dummy[k]= 1;
+ modell[k].maintype= FTYPE;
+ modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */
+ ncovf++; /* Fixed variables without age */
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 0;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else{
- printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
- fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+ printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+ fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
} /* end k1 */
}else{
printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
for(k1=1; k1<= cptcovt;k1++){
for(k2=1; k2 <k1;k2++){
if((Typevar[k1]==Typevar[k2]) && (Fixed[Tvar[k1]]==Fixed[Tvar[k2]]) && (Dummy[Tvar[k1]]==Dummy[Tvar[k2]] )){
- if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
- if(Tvar[k1]==Tvar[k2]){
- printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
- fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
- return(1);
- }
- }else if (Typevar[k1] ==2){
- k3=Tposprod[k1];
- k4=Tposprod[k2];
- if( ((Tvard[k3][1]== Tvard[k4][1])&&(Tvard[k3][2]== Tvard[k4][2])) || ((Tvard[k3][1]== Tvard[k4][2])&&(Tvard[k3][2]== Tvard[k4][1])) ){
- printf("Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
- fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
- return(1);
- }
- }
+ if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
+ if(Tvar[k1]==Tvar[k2]){
+ printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+ fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+ return(1);
+ }
+ }else if (Typevar[k1] ==2){
+ k3=Tposprod[k1];
+ k4=Tposprod[k2];
+ if( ((Tvard[k3][1]== Tvard[k4][1])&&(Tvard[k3][2]== Tvard[k4][2])) || ((Tvard[k3][1]== Tvard[k4][2])&&(Tvard[k3][2]== Tvard[k4][1])) ){
+ printf("Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+ fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+ return(1);
+ }
+ }
}
}
}
printf("ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn);
fprintf(ficlog,"ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn);
- printf("ncovf=%d, ncovv=%d, ncova=%d\n",ncovf,ncovv,ncova);
- fprintf(ficlog,"ncovf=%d, ncovv=%d, ncova=%d\n",ncovf,ncovv,ncova);
+ printf("ncovf=%d, ncovv=%d, ncova=%d, nsd=%d, nsq=%d\n",ncovf,ncovv,ncova,nsd,nsq);
+ fprintf(ficlog,"ncovf=%d, ncovv=%d, ncova=%d, nsd=%d, nsq=%d\n",ncovf,ncovv,ncova,nsd, nsq);
return (0); /* with covar[new additional covariate if product] and Tage if age */
/*endread:*/
printf("Exiting decodemodel: ");
agelim=agemaxpar;
/* i1=pow(2,ncoveff); */
- i1=pow(2,cptcoveff); /* Number of dummy covariates */
+ i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
if (cptcovn < 1){i1=1;}
for(k=1; k<=i1;k++){
char line[MAXLINE];
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
- char model[MAXLINE], modeltemp[MAXLINE];
+ char modeltemp[MAXLINE];
char resultline[MAXLINE];
char pathr[MAXLINE], pathimach[MAXLINE];
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
for(i=1; i <=nlstate; i++){
- j=0;
+ j=0;
for(jj=1; jj <=nlstate+ndeath; jj++){
- if(jj==i) continue;
- j++;
- fscanf(ficpar,"%1d%1d",&i1,&j1);
- if ((i1 != i) || (j1 != jj)){
- printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
+ if(jj==i) continue;
+ j++;
+ fscanf(ficpar,"%1d%1d",&i1,&j1);
+ if ((i1 != i) || (j1 != jj)){
+ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
It might be a problem of design; if ncovcol and the model are correct\n \
run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
- exit(1);
- }
- fprintf(ficparo,"%1d%1d",i1,j1);
- if(mle==1)
- printf("%1d%1d",i,jj);
- fprintf(ficlog,"%1d%1d",i,jj);
- for(k=1; k<=ncovmodel;k++){
- fscanf(ficpar," %lf",¶m[i][j][k]);
- if(mle==1){
- printf(" %lf",param[i][j][k]);
- fprintf(ficlog," %lf",param[i][j][k]);
- }
- else
- fprintf(ficlog," %lf",param[i][j][k]);
- fprintf(ficparo," %lf",param[i][j][k]);
- }
- fscanf(ficpar,"\n");
- numlinepar++;
- if(mle==1)
- printf("\n");
- fprintf(ficlog,"\n");
- fprintf(ficparo,"\n");
+ exit(1);
+ }
+ fprintf(ficparo,"%1d%1d",i1,j1);
+ if(mle==1)
+ printf("%1d%1d",i,jj);
+ fprintf(ficlog,"%1d%1d",i,jj);
+ for(k=1; k<=ncovmodel;k++){
+ fscanf(ficpar," %lf",¶m[i][j][k]);
+ if(mle==1){
+ printf(" %lf",param[i][j][k]);
+ fprintf(ficlog," %lf",param[i][j][k]);
+ }
+ else
+ fprintf(ficlog," %lf",param[i][j][k]);
+ fprintf(ficparo," %lf",param[i][j][k]);
+ }
+ fscanf(ficpar,"\n");
+ numlinepar++;
+ if(mle==1)
+ printf("\n");
+ fprintf(ficlog,"\n");
+ fprintf(ficparo,"\n");
}
}
fflush(ficlog);
-
+
/* Reads scales values */
p=param[1][1];
for(i=1; i <=nlstate; i++){
for(j=1; j <=nlstate+ndeath-1; j++){
- fscanf(ficpar,"%1d%1d",&i1,&j1);
- if ( (i1-i) * (j1-j) != 0){
- printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
- exit(1);
- }
- printf("%1d%1d",i,j);
- fprintf(ficparo,"%1d%1d",i1,j1);
- fprintf(ficlog,"%1d%1d",i1,j1);
- for(k=1; k<=ncovmodel;k++){
- fscanf(ficpar,"%le",&delti3[i][j][k]);
- printf(" %le",delti3[i][j][k]);
- fprintf(ficparo," %le",delti3[i][j][k]);
- fprintf(ficlog," %le",delti3[i][j][k]);
- }
- fscanf(ficpar,"\n");
- numlinepar++;
- printf("\n");
- fprintf(ficparo,"\n");
- fprintf(ficlog,"\n");
+ fscanf(ficpar,"%1d%1d",&i1,&j1);
+ if ( (i1-i) * (j1-j) != 0){
+ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
+ exit(1);
+ }
+ printf("%1d%1d",i,j);
+ fprintf(ficparo,"%1d%1d",i1,j1);
+ fprintf(ficlog,"%1d%1d",i1,j1);
+ for(k=1; k<=ncovmodel;k++){
+ fscanf(ficpar,"%le",&delti3[i][j][k]);
+ printf(" %le",delti3[i][j][k]);
+ fprintf(ficparo," %le",delti3[i][j][k]);
+ fprintf(ficlog," %le",delti3[i][j][k]);
+ }
+ fscanf(ficpar,"\n");
+ numlinepar++;
+ printf("\n");
+ fprintf(ficparo,"\n");
+ fprintf(ficlog,"\n");
}
}
fflush(ficlog);
-
+
/* Reads covariance matrix */
delti=delti3[1][1];
agedc=vector(1,n);
cod=ivector(1,n);
for(i=1;i<=n;i++){
- num[i]=0;
- moisnais[i]=0;
- annais[i]=0;
- moisdc[i]=0;
- andc[i]=0;
- agedc[i]=0;
- cod[i]=0;
- weight[i]=1.0; /* Equal weights, 1 by default */
- }
+ num[i]=0;
+ moisnais[i]=0;
+ annais[i]=0;
+ moisdc[i]=0;
+ andc[i]=0;
+ agedc[i]=0;
+ cod[i]=0;
+ weight[i]=1.0; /* Equal weights, 1 by default */
+ }
mint=matrix(1,maxwav,1,n);
anint=matrix(1,maxwav,1,n);
s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */
goto end;
/* Calculation of the number of parameters from char model */
- /* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4
+ /* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4
k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tag[cptcovage=1]=4
k=3 V4 Tvar[k=3]= 4 (from V4)
k=2 V1 Tvar[k=2]= 1 (from V1)
k=1 Tvar[1]=2 (from V2)
- */
-
- Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
+ */
+
+ Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
+ TvarsDind=ivector(1,NCOVMAX); /* */
+ TvarsD=ivector(1,NCOVMAX); /* */
+ TvarsQind=ivector(1,NCOVMAX); /* */
+ TvarsQ=ivector(1,NCOVMAX); /* */
TvarF=ivector(1,NCOVMAX); /* */
TvarFind=ivector(1,NCOVMAX); /* */
TvarV=ivector(1,NCOVMAX); /* */
free_ivector(Fixed,-1,NCOVMAX);
free_ivector(Typevar,-1,NCOVMAX);
free_ivector(Tvar,1,NCOVMAX);
+ free_ivector(TvarsQ,1,NCOVMAX);
+ free_ivector(TvarsQind,1,NCOVMAX);
+ free_ivector(TvarsD,1,NCOVMAX);
+ free_ivector(TvarsDind,1,NCOVMAX);
free_ivector(TvarFD,1,NCOVMAX);
free_ivector(TvarFDind,1,NCOVMAX);
free_ivector(TvarF,1,NCOVMAX);
git://henry.ined.fr / .git / commitdiff